Coverage for src/galsbi/ucat/plugins/apply_shear.py: 99%

1# Copyright (C) 2019 ETH Zurich, Institute for Particle Physics and Astrophysics 

2 

3""" 

4Created on Mar 27, 2019 

5author: Joerg Herbel 

6""" 

7 

8import h5py 

9import healpy as hp 

10import numpy as np 

11from cosmic_toolbox import logger 

12from ivy.plugin.base_plugin import BasePlugin 

13from ufig import coordinate_util, io_util 

14 

15from galsbi.ucat import lensing_util 

16 

17LOGGER = logger.get_logger(__file__) 

18 

19 

20def evaluate_healpix_shear_map(path, ra, dec): 

21 """ 

22 Reads in a healpix map and evaluates it at given positions. 

23 

24 :param path: path where map is stored, assumes that file contains three maps (kappa, 

25 gamma1, gamma2) 

26 :param ra: right ascension where map is evaluated 

27 :param dec: declinations where map is evaluated 

28 :return: kappa, gamma1 and gamma2 evaluated at input positions 

29 """ 

30 

31 # read map 

32 map_kappa, map_gamma_1, map_gamma_2 = hp.read_map(path, field=(0, 1, 2)) 

33 

34 # get pixel indices 

35 pixel_ind = coordinate_util.radec2pix(ra, dec, hp.npix2nside(map_kappa.size)) 

36 

37 # read out pixels 

38 kappa = map_kappa[pixel_ind] 

39 gamma_1 = map_gamma_1[pixel_ind] 

40 gamma_2 = map_gamma_2[pixel_ind] 

41 

42 return kappa, gamma_1, gamma_2 

43 

44 

45def linear_interpolation(x0, y0, x1, y1, x): 

46 """ 

47 Vectorized linear interpolation between two data points. 

48 

49 :param x0: x-coordinates of first data points 

50 :param y0: y-coordinates of first data points 

51 :param x1: x-coordinates of second data points 

52 :param y1: y-coordinates of second data points 

53 :param x: positions where interpolation is evaluated 

54 :return: interpolated values 

55 """ 

56 y = y0 + (x - x0) * (y1 - y0) / (x1 - x0) 

57 return y 

58 

59 

60def interpolate_maps_fast(z_maps, maps, z, pixel_ind, ind_intervals): 

61 """ 

62 Linearly interpolate to input redshifts between healpix maps at given redshifts. 

63 

64 :param z_maps: redshifts of input maps, assumed to be ordered 

65 :param maps: input maps corresponding to z_maps, type: hdf5-dataset or numpy-array 

66 :param z: redshifts to which maps will be interpolated, assumed to be sorted 

67 :param pixel_ind: indices indicating pixels of input maps to interpolate, same size 

68 as z 

69 :param ind_intervals: indices splitting up z into chunks that lie between the maps 

70 :return: interpolated values 

71 """ 

72 

73 pixel_unique, indices_inv = np.unique(pixel_ind, return_inverse=True) 

74 pixel_min, pixel_max = pixel_unique.min(), pixel_unique.max() 

75 maps_slab = maps[:, pixel_min : pixel_max + 1] # hdf read-in with hyperslabs 

76 maps_unique = maps_slab[:, pixel_unique - pixel_unique[0]] 

77 maps_ind = maps_unique[:, indices_inv] 

78 

79 intpol_values = np.empty_like(z) 

80 

81 # galaxies at redshifts lower than the first map --> interpolate between zero shear 

82 # at z=0 and first map 

83 intpol_values[: ind_intervals[0]] = linear_interpolation( 

84 0, 

85 np.zeros(ind_intervals[0]), 

86 z_maps[0], 

87 maps_ind[0, : ind_intervals[0]], 

88 z[: ind_intervals[0]], 

89 ) 

90 

91 # galaxies in between two maps 

92 for i_map in range(len(ind_intervals) - 1): 

93 ind_low = ind_intervals[i_map] 

94 ind_up = ind_intervals[i_map + 1] 

95 intpol_values[ind_low:ind_up] = linear_interpolation( 

96 z_maps[i_map], 

97 maps_ind[i_map, ind_low:ind_up], 

98 z_maps[i_map + 1], 

99 maps_ind[i_map + 1, ind_low:ind_up], 

100 z[ind_low:ind_up], 

101 ) 

102 

103 # galaxies at redshifts higher than the last map --> use last map to set shear 

104 # values 

105 intpol_values[ind_intervals[-1] :] = maps_ind[-1, ind_intervals[-1] :] 

106 

107 return intpol_values 

108 

109 

110def interpolate_maps(z_maps, maps, z, pixel_ind, ind_intervals): 

111 """ 

112 Linearly interpolate to input redshifts between healpix maps at given redshifts. 

113 

114 :param z_maps: redshifts of input maps, assumed to be ordered 

115 :param maps: input maps corresponding to z_maps, type: hdf5-dataset or numpy-array 

116 :param z: redshifts to which maps will be interpolated, assumed to be sorted 

117 :param pixel_ind: indices indicating pixels of input maps to interpolate, same size 

118 as z 

119 :param ind_intervals: indices splitting up z into chunks that lie between the maps 

120 :return: interpolated values 

121 """ 

122 

123 intpol_values = np.empty_like(z) 

124 map_low = np.empty(maps.shape[1], dtype=maps.dtype) 

125 map_up = maps[0][...] 

126 

127 # galaxies at redshifts lower than the first map --> interpolate between zero shear 

128 # at z=0 and first map 

129 intpol_values[: ind_intervals[0]] = linear_interpolation( 

130 0, 

131 np.zeros(ind_intervals[0]), 

132 z_maps[0], 

133 map_up[pixel_ind[: ind_intervals[0]]], 

134 z[: ind_intervals[0]], 

135 ) 

136 

137 # galaxies in between two maps 

138 for i_map in range(len(ind_intervals) - 1): 

139 ind_low = ind_intervals[i_map] 

140 ind_up = ind_intervals[i_map + 1] 

141 map_low[:] = map_up 

142 map_up[:] = maps[i_map + 1][...] 

143 

144 intpol_values[ind_low:ind_up] = linear_interpolation( 

145 z_maps[i_map], 

146 map_low[pixel_ind[ind_low:ind_up]], 

147 z_maps[i_map + 1], 

148 map_up[pixel_ind[ind_low:ind_up]], 

149 z[ind_low:ind_up], 

150 ) 

151 

152 # galaxies at redshifts higher than the last map --> use last map to set shear 

153 # values 

154 intpol_values[ind_intervals[-1] :] = map_up[pixel_ind[ind_intervals[-1] :]] 

155 

156 return intpol_values 

157 

158 

159def evaluate_hdf5_shear_maps(path, ra, dec, z): 

160 """ 

161 Evaluate hdf5 shear maps given at multiple redshifts for given redshifts and angular 

162 positions. 

163 

164 :param path: path to hdf5 file containing shear maps; assumes that this file 

165 contains four datasets: 

166 - z: redshifts of maps 

167 - kappa: kappa-maps 

168 - gamma1: gamma1-maps 

169 - gamma2: gamma2-maps 

170 :param ra: right ascensions where maps will be evaluated 

171 :param dec: declinations where maps will be evaluated 

172 :param z: redshifts to which maps will be interpolated 

173 :return: kappa, gamma1 and gamma2 

174 """ 

175 

176 # sort by redshift 

177 ind_sort = np.argsort(z) 

178 ra = ra[ind_sort] 

179 dec = dec[ind_sort] 

180 z = z[ind_sort] 

181 

182 with h5py.File(path, mode="r") as fh5: 

183 # get redshifts of maps, assumed to be sorted 

184 # print('{} reading shear map file {}'.format(str(datetime.now()), path)) 

185 LOGGER.info(f"reading shear map file {path}") 

186 z_maps = fh5["z"][...] 

187 

188 # compute pixel indices 

189 pixel_ind = coordinate_util.radec2pix( 

190 ra, dec, hp.npix2nside(fh5["kappa"].shape[1]) 

191 ) 

192 

193 # sort galaxies into redshift intervals between the maps 

194 ind_intervals = np.searchsorted(z, z_maps) 

195 

196 # read out maps and interpolate 

197 LOGGER.info(f"reading maps and interpolating for {len(z)} objects") 

198 kappa = interpolate_maps_fast(z_maps, fh5["kappa"], z, pixel_ind, ind_intervals) 

199 gamma1 = interpolate_maps_fast( 

200 z_maps, fh5["gamma1"], z, pixel_ind, ind_intervals 

201 ) 

202 gamma2 = interpolate_maps_fast( 

203 z_maps, fh5["gamma2"], z, pixel_ind, ind_intervals 

204 ) 

205 

206 # undo sorting by redshift 

207 ind_unsort = np.argsort(ind_sort) 

208 kappa[:] = kappa[ind_unsort] 

209 gamma1[:] = gamma1[ind_unsort] 

210 gamma2[:] = gamma2[ind_unsort] 

211 

212 return kappa, gamma1, gamma2 

213 

214 

215class Plugin(BasePlugin): 

216 """ 

217 Apply a potentially redshift-dependent shear from input shear maps specified by 

218 ctx.parameters.path_shear_map. Only first-order effects due to gamma are applied, 

219 kappa is completely ignored. 

220 There are three options: 

221 

222 1. The path is None, which will result in zero shear. 

223 2. The path is a file readable by healpy. The file is then assumed to contain 3 

224 healpix maps (kappa, gamma1, gamma2). 

225 3. The path is and hdf5-file containing multiple, kappa-, gamma1- and gamma2-maps at 

226 given redshifts. The shear values are computed by interpolating linearly between 

227 the maps. 

228 """ 

229 

230 def __call__(self): 

231 par = self.ctx.parameters 

232 

233 # first check if path of shear map is set to None --> zero shear 

234 if par.path_shear_map is None: 

235 self.ctx.galaxies.kappa = np.zeros(self.ctx.numgalaxies, dtype=float) 

236 self.ctx.galaxies.gamma1 = np.zeros_like(self.ctx.galaxies.kappa) 

237 self.ctx.galaxies.gamma2 = np.zeros_like(self.ctx.galaxies.kappa) 

238 

239 else: 

240 ra, dec = coordinate_util.xy2radec( 

241 coordinate_util.wcs_from_parameters(par), 

242 self.ctx.galaxies.x, 

243 self.ctx.galaxies.y, 

244 ) 

245 

246 path_shear_map = io_util.get_abs_path( 

247 par.path_shear_map, root_path=par.maps_remote_dir 

248 ) 

249 

250 if h5py.is_hdf5(path_shear_map): 

251 LOGGER.info( 

252 "Shear map file is in hdf5-format, assuming multiple" 

253 " shear maps at different redshifts" 

254 ) 

255 ( 

256 self.ctx.galaxies.kappa, 

257 self.ctx.galaxies.gamma1, 

258 self.ctx.galaxies.gamma2, 

259 ) = evaluate_hdf5_shear_maps( 

260 path=path_shear_map, ra=ra, dec=dec, z=self.ctx.galaxies.z 

261 ) 

262 else: 

263 LOGGER.info( 

264 "Shear map file is not in hdf5-format, assuming a single" 

265 " shear map stored in fits-format" 

266 ) 

267 ( 

268 self.ctx.galaxies.kappa, 

269 self.ctx.galaxies.gamma1, 

270 self.ctx.galaxies.gamma2, 

271 ) = evaluate_healpix_shear_map(path=path_shear_map, ra=ra, dec=dec) 

272 

273 # apply sign to gamma1 

274 self.ctx.galaxies.gamma1 *= par.gamma1_sign 

275 

276 # scale size 

277 self.ctx.galaxies.r50 = lensing_util.calculate_size_magnification( 

278 r=self.ctx.galaxies.int_r50, kappa=self.ctx.galaxies.kappa 

279 ) 

280 # convert to arcsec 

281 self.ctx.galaxies.r50_arcsec = self.ctx.galaxies.r50 * par.pixscale 

282 

283 # shear ellipticities 

284 LOGGER.info("applying shear and magnification") 

285 ( 

286 self.ctx.galaxies.e1, 

287 self.ctx.galaxies.e2, 

288 ) = lensing_util.apply_shear_to_ellipticities( 

289 self.ctx.galaxies.int_e1, 

290 self.ctx.galaxies.int_e2, 

291 self.ctx.galaxies.kappa, 

292 self.ctx.galaxies.gamma1, 

293 self.ctx.galaxies.gamma2, 

294 ) 

295 

296 # magnify flux 

297 

298 # magnification in terms of flux, see Bartelmann & Schneider 2001, 

299 # https://arxiv.org/pdf/astro-ph/9912508.pdf, eq. (3.14) 

300 magnification = lensing_util.calculate_flux_magnification( 

301 self.ctx.galaxies.kappa, self.ctx.galaxies.gamma1, self.ctx.galaxies.gamma2 

302 ) 

303 

304 # in terms of magnitudes 

305 magnification[:] = 2.5 * np.log10(magnification) 

306 

307 for mag in self.ctx.galaxies.magnitude_dict.values(): 

308 mag -= magnification 

309 

310 # add new columns to list 

311 self.ctx.galaxies.columns.extend( 

312 ["gamma1", "gamma2", "kappa", "e1", "e2", "r50", "r50_arcsec"] 

313 ) 

314 

315 def __str__(self): 

316 return "apply shear"