Coverage for src/galsbi/ucat/plugins/apply_shear.py: 99%
89 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-13 03:24 +0000
1# Copyright (C) 2019 ETH Zurich, Institute for Particle Physics and Astrophysics
3"""
4Created on Mar 27, 2019
5author: Joerg Herbel
6"""
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
15from galsbi.ucat import lensing_util
17LOGGER = logger.get_logger(__file__)
20def evaluate_healpix_shear_map(path, ra, dec):
21 """
22 Reads in a healpix map and evaluates it at given positions.
23 :param path: path where map is stored, assumes that file contains three maps (kappa,
24 gamma1, gamma2)
25 :param ra: right ascension where map is evaluated
26 :param dec: declinations where map is evaluated
27 :return: kappa, gamma1 and gamma2 evaluated at input positions
28 """
30 # read map
31 map_kappa, map_gamma_1, map_gamma_2 = hp.read_map(path, field=(0, 1, 2))
33 # get pixel indices
34 pixel_ind = coordinate_util.radec2pix(ra, dec, hp.npix2nside(map_kappa.size))
36 # read out pixels
37 kappa = map_kappa[pixel_ind]
38 gamma_1 = map_gamma_1[pixel_ind]
39 gamma_2 = map_gamma_2[pixel_ind]
41 return kappa, gamma_1, gamma_2
44def linear_interpolation(x0, y0, x1, y1, x):
45 """
46 Vectorized linear interpolation between two data points.
47 :param x0: x-coordinates of first data points
48 :param y0: y-coordinates of first data points
49 :param x1: x-coordinates of second data points
50 :param y1: y-coordinates of second data points
51 :param x: positions where interpolation is evaluated
52 :return: interpolated values
53 """
54 y = y0 + (x - x0) * (y1 - y0) / (x1 - x0)
55 return y
58def interpolate_maps_fast(z_maps, maps, z, pixel_ind, ind_intervals):
59 """
60 Linearly interpolate to input redshifts between healpix maps at given redshifts.
61 :param z_maps: redshifts of input maps, assumed to be ordered
62 :param maps: input maps corresponding to z_maps, type: hdf5-dataset or numpy-array
63 :param z: redshifts to which maps will be interpolated, assumed to be sorted
64 :param pixel_ind: indices indicating pixels of input maps to interpolate, same size
65 as z
66 :param ind_intervals: indices splitting up z into chunks that lie between the maps
67 :return: interpolated values
68 """
70 pixel_unique, indices_inv = np.unique(pixel_ind, return_inverse=True)
71 pixel_min, pixel_max = pixel_unique.min(), pixel_unique.max()
72 maps_slab = maps[:, pixel_min : pixel_max + 1] # hdf read-in with hyperslabs
73 maps_unique = maps_slab[:, pixel_unique - pixel_unique[0]]
74 maps_ind = maps_unique[:, indices_inv]
76 intpol_values = np.empty_like(z)
78 # galaxies at redshifts lower than the first map --> interpolate between zero shear
79 # at z=0 and first map
80 intpol_values[: ind_intervals[0]] = linear_interpolation(
81 0,
82 np.zeros(ind_intervals[0]),
83 z_maps[0],
84 maps_ind[0, : ind_intervals[0]],
85 z[: ind_intervals[0]],
86 )
88 # galaxies in between two maps
89 for i_map in range(len(ind_intervals) - 1):
90 ind_low = ind_intervals[i_map]
91 ind_up = ind_intervals[i_map + 1]
92 intpol_values[ind_low:ind_up] = linear_interpolation(
93 z_maps[i_map],
94 maps_ind[i_map, ind_low:ind_up],
95 z_maps[i_map + 1],
96 maps_ind[i_map + 1, ind_low:ind_up],
97 z[ind_low:ind_up],
98 )
100 # galaxies at redshifts higher than the last map --> use last map to set shear
101 # values
102 intpol_values[ind_intervals[-1] :] = maps_ind[-1, ind_intervals[-1] :]
104 return intpol_values
107def interpolate_maps(z_maps, maps, z, pixel_ind, ind_intervals):
108 """
109 Linearly interpolate to input redshifts between healpix maps at given redshifts.
110 :param z_maps: redshifts of input maps, assumed to be ordered
111 :param maps: input maps corresponding to z_maps, type: hdf5-dataset or numpy-array
112 :param z: redshifts to which maps will be interpolated, assumed to be sorted
113 :param pixel_ind: indices indicating pixels of input maps to interpolate, same size
114 as z
115 :param ind_intervals: indices splitting up z into chunks that lie between the maps
116 :return: interpolated values
117 """
119 intpol_values = np.empty_like(z)
120 map_low = np.empty(maps.shape[1], dtype=maps.dtype)
121 map_up = maps[0][...]
123 # galaxies at redshifts lower than the first map --> interpolate between zero shear
124 # at z=0 and first map
125 intpol_values[: ind_intervals[0]] = linear_interpolation(
126 0,
127 np.zeros(ind_intervals[0]),
128 z_maps[0],
129 map_up[pixel_ind[: ind_intervals[0]]],
130 z[: ind_intervals[0]],
131 )
133 # galaxies in between two maps
134 for i_map in range(len(ind_intervals) - 1):
135 ind_low = ind_intervals[i_map]
136 ind_up = ind_intervals[i_map + 1]
137 map_low[:] = map_up
138 map_up[:] = maps[i_map + 1][...]
140 intpol_values[ind_low:ind_up] = linear_interpolation(
141 z_maps[i_map],
142 map_low[pixel_ind[ind_low:ind_up]],
143 z_maps[i_map + 1],
144 map_up[pixel_ind[ind_low:ind_up]],
145 z[ind_low:ind_up],
146 )
148 # galaxies at redshifts higher than the last map --> use last map to set shear
149 # values
150 intpol_values[ind_intervals[-1] :] = map_up[pixel_ind[ind_intervals[-1] :]]
152 return intpol_values
155def evaluate_hdf5_shear_maps(path, ra, dec, z):
156 """
157 Evaluate hdf5 shear maps given at multiple redshifts for given redshifts and angular
158 positions.
159 :param path: path to hdf5 file containing shear maps; assumes that this file
160 contains four datasets:
161 - z: redshifts of maps
162 - kappa: kappa-maps
163 - gamma1: gamma1-maps
164 - gamma2: gamma2-maps
165 :param ra: right ascensions where maps will be evaluated
166 :param dec: declinations where maps will be evaluated
167 :param z: redshifts to which maps will be interpolated
168 :return: kappa, gamma1 and gamma2
169 """
171 # sort by redshift
172 ind_sort = np.argsort(z)
173 ra = ra[ind_sort]
174 dec = dec[ind_sort]
175 z = z[ind_sort]
177 with h5py.File(path, mode="r") as fh5:
178 # get redshifts of maps, assumed to be sorted
179 # print('{} reading shear map file {}'.format(str(datetime.now()), path))
180 LOGGER.info(f"reading shear map file {path}")
181 z_maps = fh5["z"][...]
183 # compute pixel indices
184 pixel_ind = coordinate_util.radec2pix(
185 ra, dec, hp.npix2nside(fh5["kappa"].shape[1])
186 )
188 # sort galaxies into redshift intervals between the maps
189 ind_intervals = np.searchsorted(z, z_maps)
191 # read out maps and interpolate
192 LOGGER.info(f"reading maps and interpolating for {len(z)} objects")
193 kappa = interpolate_maps_fast(z_maps, fh5["kappa"], z, pixel_ind, ind_intervals)
194 gamma1 = interpolate_maps_fast(
195 z_maps, fh5["gamma1"], z, pixel_ind, ind_intervals
196 )
197 gamma2 = interpolate_maps_fast(
198 z_maps, fh5["gamma2"], z, pixel_ind, ind_intervals
199 )
201 # undo sorting by redshift
202 ind_unsort = np.argsort(ind_sort)
203 kappa[:] = kappa[ind_unsort]
204 gamma1[:] = gamma1[ind_unsort]
205 gamma2[:] = gamma2[ind_unsort]
207 return kappa, gamma1, gamma2
210class Plugin(BasePlugin):
211 """
212 Apply a potentially redshift-dependent shear from input shear maps specified by
213 ctx.parameters.path_shear_map. Only
214 first-order effects due to gamma are applied, kappa is completely ignored.
215 There are three options:
216 1) The path is None, which will result in zero shear.
217 2) The path is a file readable by healpy. The file is then assumed to contain 3
218 healpix maps (kappa, gamma1, gamma2).
219 3) The path is and hdf5-file containing multiple, kappa-, gamma1- and gamma2-maps at
220 given redshifts. The shear values are computed by interpolating linearly between
221 the maps.
222 """
224 def __call__(self):
225 par = self.ctx.parameters
227 # first check if path of shear map is set to None --> zero shear
228 if par.path_shear_map is None:
229 self.ctx.galaxies.kappa = np.zeros(self.ctx.numgalaxies, dtype=float)
230 self.ctx.galaxies.gamma1 = np.zeros_like(self.ctx.galaxies.kappa)
231 self.ctx.galaxies.gamma2 = np.zeros_like(self.ctx.galaxies.kappa)
233 else:
234 ra, dec = coordinate_util.xy2radec(
235 coordinate_util.wcs_from_parameters(par),
236 self.ctx.galaxies.x,
237 self.ctx.galaxies.y,
238 )
240 path_shear_map = io_util.get_abs_path(
241 par.path_shear_map, root_path=par.maps_remote_dir
242 )
244 if h5py.is_hdf5(path_shear_map):
245 LOGGER.info(
246 "Shear map file is in hdf5-format, assuming multiple"
247 " shear maps at different redshifts"
248 )
249 (
250 self.ctx.galaxies.kappa,
251 self.ctx.galaxies.gamma1,
252 self.ctx.galaxies.gamma2,
253 ) = evaluate_hdf5_shear_maps(
254 path=path_shear_map, ra=ra, dec=dec, z=self.ctx.galaxies.z
255 )
256 else:
257 LOGGER.info(
258 "Shear map file is not in hdf5-format, assuming a single"
259 " shear map stored in fits-format"
260 )
261 (
262 self.ctx.galaxies.kappa,
263 self.ctx.galaxies.gamma1,
264 self.ctx.galaxies.gamma2,
265 ) = evaluate_healpix_shear_map(path=path_shear_map, ra=ra, dec=dec)
267 # apply sign to gamma1
268 self.ctx.galaxies.gamma1 *= par.gamma1_sign
270 # scale size
271 self.ctx.galaxies.r50 = lensing_util.calculate_size_magnification(
272 r=self.ctx.galaxies.int_r50, kappa=self.ctx.galaxies.kappa
273 )
275 # shear ellipticities
276 LOGGER.info("applying shear and magnification")
277 (
278 self.ctx.galaxies.e1,
279 self.ctx.galaxies.e2,
280 ) = lensing_util.apply_shear_to_ellipticities(
281 self.ctx.galaxies.int_e1,
282 self.ctx.galaxies.int_e2,
283 self.ctx.galaxies.kappa,
284 self.ctx.galaxies.gamma1,
285 self.ctx.galaxies.gamma2,
286 )
288 # magnify flux
290 # magnification in terms of flux, see Bartelmann & Schneider 2001,
291 # https://arxiv.org/pdf/astro-ph/9912508.pdf, eq. (3.14)
292 magnification = lensing_util.calculate_flux_magnification(
293 self.ctx.galaxies.kappa, self.ctx.galaxies.gamma1, self.ctx.galaxies.gamma2
294 )
296 # in terms of magnitudes
297 magnification[:] = 2.5 * np.log10(magnification)
299 for mag in self.ctx.galaxies.magnitude_dict.values():
300 mag -= magnification
302 # add new columns to list
303 self.ctx.galaxies.columns.extend(
304 ["gamma1", "gamma2", "kappa", "e1", "e2", "r50"]
305 )
307 def __str__(self):
308 return "apply shear"