Coverage for metadetect / lsst / metacal_exposures.py: 12%

1""" 

2Code to do metacal with lsst exposures 

3""" 

4import numpy as np 

5from ngmix.metacal.metacal import _get_gauss_target_psf 

6import galsim 

7import lsst.afw.image as afw_image 

8from .util import ( 

9 get_integer_center, get_jacobian, get_stack_kernel_psf, get_mbexp, 

10) 

11 

12DEFAULT_TYPES = ['noshear', '1p', '1m'] 

13INTERP = 'lanczos15' 

14STEP = 0.01 

15 

16 

17def get_metacal_mbexps_fixnoise(mbexp, noise_mbexp, types=None): 

18 """ 

19 Get metacal MultibandExposures with fixed noise 

20 

21 Parameters 

22 ---------- 

23 mbexp: lsst.afw.image.MultibandExposure 

24 The exposure data 

25 noise_mbexp: lsst.afw.image.MultibandExposure 

26 The exposure data with pure noise 

27 types: list, optional 

28 The metacal types, e.g. ('noshear', '1p', '1m') 

29 

30 Returns 

31 ------- 

32 mdict, noise_mdict 

33 dicts keyed by type, holding exposures 

34 """ 

35 

36 mdict = get_metacal_mbexps(mbexp=mbexp, types=types) 

37 noise_mdict = get_metacal_mbexps(mbexp=noise_mbexp, types=types, rot=True) 

38 for shear_type in mdict: 

39 for exp, nexp in zip(mdict[shear_type], noise_mdict[shear_type]): 

40 exp.image.array[:, :] += nexp.image.array[:, :] 

41 exp.variance.array[:, :] *= 2 

42 

43 return mdict, noise_mdict 

44 

45 

46def get_metacal_mbexps(mbexp, types=None, rot=False): 

47 """ 

48 Get metacal MultibandExposures 

49 

50 Parameters 

51 ---------- 

52 mbexp: lsst.afw.image.MultibandExposure 

53 The exposure data 

54 types: list, optional 

55 The metacal types, e.g. ('noshear', '1p', '1m') 

56 rot: bool, optional 

57 If set to True, rotate before shearing, then rotate back. 

58 

59 Returns 

60 ------- 

61 mdict 

62 dict keyed by type, holding exposures 

63 """ 

64 

65 if types is None: 

66 types = DEFAULT_TYPES 

67 

68 mdict_with_explists = {} 

69 for shear_type in types: 

70 mdict_with_explists[shear_type] = [] 

71 

72 for iband, band in enumerate(mbexp.bands): 

73 exp = mbexp[band] 

74 

75 this_mdict = get_metacal_exps(exp, types=types, rot=rot) 

76 

77 for shear_type in this_mdict: 

78 exp = this_mdict[shear_type] 

79 mdict_with_explists[shear_type].append(exp) 

80 

81 # now build the mbexp 

82 mdict = {} 

83 for shear_type in types: 

84 # this properly copies over the wcs and filter label 

85 mdict[shear_type] = get_mbexp(mdict_with_explists[shear_type]) 

86 

87 return mdict 

88 

89 

90def get_metacal_exps_fixnoise(exp, noise_exp, types=None): 

91 """ 

92 Get metacal exposures with fixed noise 

93 

94 Parameters 

95 ---------- 

96 exp: lsst.afw.image.Exposure 

97 The exposure data 

98 noise_exp: lsst.afw.image.Exposure 

99 The exposure data with pure noise 

100 types: list, optional 

101 The metacal types, e.g. ('noshear', '1p', '1m') 

102 

103 Returns 

104 ------- 

105 dict keyed by type, holding exposures 

106 """ 

107 if types is None: 

108 types = DEFAULT_TYPES 

109 

110 mdict = get_metacal_exps(exp, types=types) 

111 noise_mdict = get_metacal_exps(noise_exp, types=types, rot=True) 

112 

113 for shear_type in types: 

114 exp = mdict[shear_type] 

115 nexp = noise_mdict[shear_type] 

116 

117 exp.image.array[:, :] += nexp.image.array[:, :] 

118 exp.variance.array[:, :] *= 2 

119 

120 return mdict, noise_mdict 

121 

122 

123def get_metacal_exps(exp, types=None, rot=False): 

124 """ 

125 Get metacal exposures 

126 

127 Parameters 

128 ---------- 

129 exp: lsst.afw.image.Exposure 

130 The exposure data 

131 types: list, optional 

132 The metacal types, e.g. ('noshear', '1p', '1m') 

133 rot: bool, optional 

134 If set to True, rotate before shearing, then rotate back. 

135 

136 Returns 

137 ------- 

138 dict keyed by type, holding exposures 

139 """ 

140 

141 if types is None: 

142 types = DEFAULT_TYPES 

143 

144 cen, _ = get_integer_center(exp.getWcs(), exp.getBBox(), as_double=True) 

145 

146 gwcs = get_galsim_jacobian_wcs(exp=exp, cen=cen) 

147 

148 psf_image_array = get_psf_kernel_image(exp=exp, cen=cen) 

149 psf_flux = psf_image_array.sum() 

150 

151 eimage = exp.image.array.copy() 

152 if rot: 

153 eimage = np.rot90(eimage, k=1) 

154 

155 # force double, *it matters*. There are clear artifacts 

156 # present when comparing to the regular metacal if float32 is 

157 # used 

158 image = galsim.ImageD(eimage, wcs=gwcs) 

159 psf_image = galsim.ImageD(psf_image_array, wcs=gwcs) 

160 

161 # interpolated psf image and its deconvolution 

162 psf_int = galsim.InterpolatedImage(psf_image, x_interpolant=INTERP) 

163 

164 # deconvolved galaxy image, psf+pixel removed 

165 image_int_nopsf = galsim.Convolve( 

166 galsim.InterpolatedImage(image, x_interpolant=INTERP), 

167 galsim.Deconvolve(psf_int), 

168 ) 

169 

170 gauss_psf = _get_gauss_target_psf(psf_int, flux=psf_flux) 

171 

172 dilation = 1.0 + 2.0*STEP 

173 psf_dilated = gauss_psf.dilate(dilation) 

174 

175 psf_dilated_image = psf_image.copy() 

176 psf_dilated.drawImage(image=psf_dilated_image, method='no_pixel') 

177 

178 mdict = {} 

179 for shear_type in types: 

180 mdict[shear_type] = _get_metacal_exp( 

181 exp=exp, image=image, image_int_nopsf=image_int_nopsf, 

182 psf_dilated=psf_dilated, psf_dilated_image=psf_dilated_image, 

183 shear_type=shear_type, rot=rot, 

184 ) 

185 

186 return mdict 

187 

188 

189def _get_metacal_exp( 

190 exp, image, image_int_nopsf, psf_dilated, psf_dilated_image, shear_type, 

191 rot, 

192): 

193 sexp = afw_image.ExposureD(exp, deep=True) 

194 

195 if shear_type == 'noshear': 

196 obj = image_int_nopsf 

197 else: 

198 shear = get_shear_from_type(shear_type) 

199 obj = image_int_nopsf.shear(shear) 

200 

201 shimage = image.copy() 

202 galsim.Convolve(obj, psf_dilated).drawImage( 

203 image=shimage, method='no_pixel', 

204 ) 

205 

206 simage = shimage.array 

207 if rot: 

208 simage = np.rot90(simage, k=3) 

209 

210 sexp.image.array[:, :] = simage 

211 

212 stack_psf = get_stack_kernel_psf(psf_dilated_image.array) 

213 sexp.setPsf(stack_psf) 

214 return sexp 

215 

216 

217def get_shear_from_type(shear_type): 

218 """ 

219 Convert a shear type string to a galsim Shear 

220 

221 Parameters 

222 ---------- 

223 shear_type: str 

224 One of ('1p', '1m', '2p', '1m') 

225 

226 Returns 

227 ------- 

228 galsim.Shear 

229 """ 

230 if shear_type == '1p': 

231 shear = galsim.Shear(g1=STEP) 

232 elif shear_type == '1m': 

233 shear = galsim.Shear(g1=-STEP) 

234 elif shear_type == '2p': 

235 shear = galsim.Shear(g2=STEP) 

236 elif shear_type == '2m': 

237 shear = galsim.Shear(g2=-STEP) 

238 else: 

239 raise ValueError('shear type should be noshear, 1p, 1m, 2p, 2m') 

240 

241 return shear 

242 

243 

244def get_galsim_jacobian_wcs(exp, cen): 

245 """ 

246 Get a galsim jacobian wcs object from the input exposure's 

247 wcs and the input location 

248 

249 Parameters 

250 ----------- 

251 exp: lsst.afw.image.Exposure* 

252 An lsst exposure object 

253 cen: lsst.afw.geom.Point2D 

254 The point in the image at which to construct the wcs 

255 

256 Returns 

257 ------- 

258 galsim.JacobianWCS 

259 """ 

260 jac = get_jacobian(exp=exp, cen=cen) 

261 return jac.get_galsim_wcs() 

262 

263 

264def get_psf_kernel_image(exp, cen): 

265 """ 

266 Reconstruct a kernel image for the input exposure at the specified 

267 location 

268 

269 Parameters 

270 ----------- 

271 exp: lsst.afw.image.Exposure* 

272 An lsst exposure object 

273 cen: lsst.afw.geom.Point2D 

274 The point in the image at which to construct the wcs 

275 

276 Returns 

277 ------- 

278 image as a numpy array 

279 """ 

280 psf_obj = exp.getPsf() 

281 return psf_obj.computeKernelImage(cen).array