Coverage for python / lsst / fgcmcal / focalPlaneProjector.py: 19%

1# This file is part of fgcmcal. 

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"""A class to project the focal plane in arbitrary rotations for fgcm. 

22 

23This file contains a class used by fgcm ... 

24""" 

25from functools import lru_cache 

26import warnings 

27import numpy as np 

28 

29import lsst.afw.image as afwImage 

30import lsst.afw.cameraGeom as afwCameraGeom 

31import lsst.geom as geom 

32from lsst.obs.base import createInitialSkyWcs 

33 

34__all__ = ['FocalPlaneProjector'] 

35 

36 

37class FocalPlaneProjector(object): 

38 """ 

39 Class to project the focal plane onto the sky. 

40 

41 Parameters 

42 ---------- 

43 camera : `lsst.afw.cameraGeom.Camera` 

44 Camera from the butler. 

45 defaultOrientation : `int` 

46 Default camera orientation in degrees. This angle is the position 

47 angle of the focal plane +Y with respect to north. 

48 useScienceDetectors : `bool`, optional 

49 Use only science detectors in projector? 

50 """ 

51 def __init__(self, camera, defaultOrientation, useScienceDetectors=False): 

52 self.camera = camera 

53 

54 # Put the reference boresight at the equator to avoid cos(dec) problems. 

55 self.boresight = geom.SpherePoint(180.0*geom.degrees, 0.0*geom.degrees) 

56 self.flipX = False 

57 self.defaultOrientation = int(defaultOrientation) % 360 

58 self.useScienceDetectors = useScienceDetectors 

59 

60 def _makeWcsDict(self, orientation): 

61 """ 

62 Make a dictionary of WCSs at the reference boresight position. 

63 

64 Parameters 

65 ---------- 

66 orientation : `int` 

67 Orientation in degrees. This angle is the position 

68 angle of the focal plane +Y with respect to north. 

69 

70 Returns 

71 ------- 

72 wcsDict : `dict` 

73 Dictionary of WCS, with the detector id as the key. 

74 """ 

75 _orientation = orientation*geom.degrees 

76 

77 visitInfo = afwImage.VisitInfo(boresightRaDec=self.boresight, 

78 boresightRotAngle=_orientation, 

79 rotType=afwImage.RotType.SKY) 

80 

81 wcsDict = {} 

82 

83 for detector in self.camera: 

84 if self.useScienceDetectors: 

85 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

86 continue 

87 

88 detectorId = detector.getId() 

89 wcsDict[detectorId] = createInitialSkyWcs(visitInfo, detector, self.flipX) 

90 

91 return wcsDict 

92 

93 def __call__(self, orientation, nstep=100, use_cache=True): 

94 """ 

95 Make a focal plane projection mapping for use with fgcm. 

96 

97 Parameters 

98 ---------- 

99 orientation : `float` or `int` 

100 Camera orientation in degrees. This angle is the position 

101 angle of the focal plane +Y with respect to north. 

102 nstep : `int` 

103 Number of steps in x/y per detector for the mapping. 

104 use_cache : `bool`, optional 

105 Use integerized cached lookup. 

106 

107 Returns 

108 ------- 

109 projectionMapping : `np.ndarray` 

110 A projection mapping object with x, y, x_size, y_size, 

111 delta_ra_cent, delta_dec_cent, delta_ra, delta_dec for 

112 each detector id. 

113 """ 

114 if not np.isfinite(orientation): 

115 warnings.warn('Encountered non-finite orientation; using default.') 

116 _orientation = self.defaultOrientation 

117 else: 

118 _orientation = orientation % 360 

119 

120 if use_cache: 

121 _orientation = int(_orientation) 

122 

123 return self._compute_cached_projection(int(_orientation), nstep=nstep) 

124 else: 

125 return self._compute_projection(_orientation, nstep=nstep) 

126 

127 @lru_cache(maxsize=360) 

128 def _compute_cached_projection(self, orientation, nstep=50): 

129 """ 

130 Compute the focal plane projection, with caching. 

131 

132 Parameters 

133 ---------- 

134 orientation : `int` 

135 Camera orientation in degrees. This angle is the position 

136 angle of the focal plane +Y with respect to north. 

137 nstep : `int` 

138 Number of steps in x/y per detector for the mapping. 

139 

140 Returns 

141 ------- 

142 projectionMapping : `np.ndarray` 

143 A projection mapping object with x, y, x_size, y_size, 

144 delta_ra_cent, delta_dec_cent, delta_ra, delta_dec for 

145 each detector id. 

146 """ 

147 return self._compute_projection(orientation, nstep=nstep) 

148 

149 def _compute_projection(self, orientation, nstep=50): 

150 """ 

151 Compute the focal plane projection. 

152 

153 Parameters 

154 ---------- 

155 orientation : `float` or `int` 

156 Camera orientation in degrees. This angle is the position 

157 angle of the focal plane +Y with respect to north. 

158 nstep : `int` 

159 Number of steps in x/y per detector for the mapping. 

160 

161 Returns 

162 ------- 

163 projectionMapping : `np.ndarray` 

164 A projection mapping object with x, y, x_size, y_size, 

165 delta_ra_cent, delta_dec_cent, delta_ra, delta_dec for 

166 each detector id. 

167 """ 

168 wcsDict = self._makeWcsDict(orientation) 

169 

170 # Need something for the max detector ... 

171 deltaMapper = np.zeros( 

172 len(wcsDict), 

173 dtype=[ 

174 ('id', 'i4'), 

175 ('x', 'f8', nstep**2), 

176 ('y', 'f8', nstep**2), 

177 ('x_size', 'i4'), 

178 ('y_size', 'i4'), 

179 ('delta_ra_cent', 'f8'), 

180 ('delta_dec_cent', 'f8'), 

181 ('delta_ra', 'f8', nstep**2), 

182 ('delta_dec', 'f8', nstep**2) 

183 ], 

184 ) 

185 

186 for detector in self.camera: 

187 if self.useScienceDetectors: 

188 if not detector.getType() == afwCameraGeom.DetectorType.SCIENCE: 

189 continue 

190 

191 detectorId = detector.getId() 

192 

193 deltaMapper['id'][detectorId] = detectorId 

194 

195 xSize = detector.getBBox().getMaxX() 

196 ySize = detector.getBBox().getMaxY() 

197 

198 xValues = np.linspace(0.0, xSize, nstep) 

199 yValues = np.linspace(0.0, ySize, nstep) 

200 

201 deltaMapper['x'][detectorId, :] = np.repeat(xValues, yValues.size) 

202 deltaMapper['y'][detectorId, :] = np.tile(yValues, xValues.size) 

203 deltaMapper['x_size'][detectorId] = xSize 

204 deltaMapper['y_size'][detectorId] = ySize 

205 

206 radec = wcsDict[detector.getId()].pixelToSkyArray(deltaMapper['x'][detectorId, :], 

207 deltaMapper['y'][detectorId, :], 

208 degrees=True) 

209 

210 deltaMapper['delta_ra'][detectorId, :] = radec[0] - self.boresight.getRa().asDegrees() 

211 deltaMapper['delta_dec'][detectorId, :] = radec[1] - self.boresight.getDec().asDegrees() 

212 

213 detCenter = wcsDict[detector.getId()].pixelToSky(detector.getCenter(afwCameraGeom.PIXELS)) 

214 deltaMapper['delta_ra_cent'][detectorId] = (detCenter.getRa() 

215 - self.boresight.getRa()).asDegrees() 

216 deltaMapper['delta_dec_cent'][detectorId] = (detCenter.getDec() 

217 - self.boresight.getDec()).asDegrees() 

218 

219 return deltaMapper