Coverage for src/ufig/plugins/draw_stars_besancon_map.py: 94%

178 statements  

« prev     ^ index     » next       coverage.py v7.6.9, created at 2024-12-12 19:08 +0000

1# Copyright (C) 2016 ETH Zurich, Institute for Astronomy 

2 

3""" 

4Created on Jul 20, 2016 

5author: Tomasz Kacprzak 

6""" 

7 

8import warnings 

9 

10import h5py 

11import healpy as hp 

12import numpy as np 

13import six 

14from cosmic_toolbox import arraytools as at 

15from ivy.plugin.base_plugin import BasePlugin 

16 

17from .. import coordinate_util, io_util, sampling_util 

18from ..sampling_util import sample_position_uniform 

19 

20if six.PY2: # pragma: no cover 

21 import cPickle as pickle 

22else: 

23 import pickle 

24 

25 

26NAME = "draw stars" 

27 

28 

29def get_interp_nearest(ra_new, dec_new, dict_besancon_info): 

30 """ 

31 

32 :param ra_new: 

33 :param dec_new: 

34 :param dict_besancon_info: 

35 :return: 

36 """ 

37 

38 ipix = coordinate_util.radec2pix( 

39 ra=ra_new, dec=dec_new, nside=dict_besancon_info["nside"] 

40 ) 

41 if dict_besancon_info["healpix_list"][ipix] is None: 41 ↛ 42line 41 didn't jump to line 42 because the condition on line 41 was never true

42 warnings.warn( 

43 "Besancon model: queried position is outside the area covered by the model;" 

44 " falling back to default", 

45 stacklevel=1, 

46 ) 

47 ipix = 10000 

48 

49 closest_cat = dict_besancon_info["healpix_list"][ipix] 

50 

51 return closest_cat 

52 

53 

54def load_besancon_map_info(ctx): 

55 par = ctx.parameters 

56 

57 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir) 

58 

59 dict_besancon_info = {} 

60 

61 with h5py.File(filepath_map, "r") as ff: 

62 dict_besancon_info["simulation_area"] = float(np.array(ff["simulation_area"])) 

63 dict_besancon_info["nside"] = int(np.array(ff["nside"])) 

64 

65 return dict_besancon_info 

66 

67 

68def load_besancon_map_pixels(ctx, list_ipix): 

69 """ 

70 Simply load the pickle with maps, from par.besancon_map_path 

71 :return pickle with the Besancon simulation map model, with fields: 

72 

73 dict_besancon_info={'simulation_area': area_store, 'nside': nside, 

74 'healpix_list': list_df, 'healpix_mask': hp_map} 

75 

76 simulation_area - area in sq. deg corresponding to the simulation catalog covers 

77 

78 nside - resolution of simulation sampling 

79 

80 healpix_list - a list of numpy record arrays, each element corresponds to a pixel on 

81 HEALPix grid with nside=nside 

82 

83 hp_map - HEALPix map showing the coverage of simulations, ring=False 

84 """ 

85 

86 par = ctx.parameters 

87 

88 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir) 

89 

90 dist_besancon_pix_cat = {} 

91 

92 with h5py.File(filepath_map, "r") as ff: 

93 for ipix in list_ipix: 

94 if len(ff[f"healpix_list/{ipix:04d}"]) == 0: 94 ↛ 95line 94 didn't jump to line 95 because the condition on line 94 was never true

95 besancon_pix_cat = None 

96 else: 

97 besancon_pix_cat = np.array(ff[f"healpix_list/{ipix:04d}"]) 

98 

99 dist_besancon_pix_cat[f"{ipix:04d}"] = besancon_pix_cat 

100 

101 return dist_besancon_pix_cat 

102 

103 

104def load_besancon_map(ctx): 

105 """ 

106 Simply load the pickle with maps, from par.besancon_map_path 

107 :return pickle with the Besancon simulation map model, with fields: 

108 

109 dict_besancon_info = {'simulation_area': area_store, 'nside': nside, 

110 'healpix_list': list_df, 'healpix_mask': hp_map} 

111 

112 simulation_area - area in sq. deg corresponding to the simulation catalog covers 

113 

114 nside - resolution of simulation sampling 

115 

116 healpix_list - a list of numpy record arrays, each element corresponds 

117 to a pixel on HEALPix grid with nside=nside 

118 

119 hp_map - HEALPix map showing the coverage of simulations, ring=False 

120 """ 

121 

122 par = ctx.parameters 

123 

124 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir) 

125 

126 if ".pkl" in par.besancon_map_path: 126 ↛ 127line 126 didn't jump to line 127 because the condition on line 126 was never true

127 with open(filepath_map) as ff: 

128 dict_besancon_info = pickle.load(ff) 

129 

130 elif ".h5" in par.besancon_map_path: 130 ↛ 148line 130 didn't jump to line 148 because the condition on line 130 was always true

131 dict_besancon_info = {} 

132 with h5py.File(filepath_map, "r") as ff: 

133 dict_besancon_info["healpix_mask"] = np.array(ff["healpix_mask"]) 

134 dict_besancon_info["simulation_area"] = float( 

135 np.array(ff["simulation_area"]) 

136 ) 

137 dict_besancon_info["nside"] = int(np.array(ff["nside"])) 

138 n_pix = len(dict_besancon_info["healpix_mask"]) 

139 dict_besancon_info["healpix_list"] = [None] * n_pix 

140 for ip in range(n_pix): 

141 if len(ff[f"healpix_list/{ip:04d}"]) == 0: 

142 dict_besancon_info["healpix_list"][ip] = None 

143 else: 

144 dict_besancon_info["healpix_list"][ip] = np.array( 

145 ff[f"healpix_list/{ip:04d}"] 

146 ) 

147 

148 return dict_besancon_info 

149 

150 

151def get_star_cat_besancon(ctx): 

152 par = ctx.parameters 

153 

154 w = coordinate_util.tile_in_skycoords( 

155 pixscale=par.pixscale, 

156 ra0=par.ra0, 

157 dec0=par.dec0, 

158 crpix_ra=par.crpix_ra, 

159 crpix_dec=par.crpix_dec, 

160 ) 

161 

162 pixels = coordinate_util.get_healpix_pixels( 

163 par.nside_sampling, w, par.size_x, par.size_y 

164 ) 

165 pixarea = hp.nside2pixarea(par.nside_sampling, degrees=True) 

166 

167 dict_besancon_info = load_besancon_map(ctx) 

168 cat_stars = get_interp_nearest( 

169 ra_new=par.ra0, dec_new=par.dec0, dict_besancon_info=dict_besancon_info 

170 ) 

171 simulation_area = dict_besancon_info["simulation_area"] 

172 

173 # Choose stars randomly according to magnitude limits (including Poisson noise), 

174 # apply limits for detection band 

175 mask = (cat_stars[par.reference_band] >= par.stars_mag_min) & ( 

176 cat_stars[par.reference_band] <= par.stars_mag_max 

177 ) 

178 

179 cat_stars_select = cat_stars[mask] 

180 

181 max_numstars = np.count_nonzero(mask) 

182 

183 # Mean number of stars within a Healpix pixel 

184 mean_numstars = np.int32(np.around(max_numstars * pixarea / simulation_area)) 

185 

186 # Loop over pixels 

187 for pixel in pixels: 

188 # Reseed random library 

189 np.random.seed(par.seed + pixel + par.star_num_seed_offset) 

190 nstars = np.random.poisson(mean_numstars) 

191 while ( 191 ↛ 194line 191 didn't jump to line 194

192 nstars > max_numstars 

193 ): # in case drawn number of stars is larger than number of stars in catalog 

194 nstars = np.random.poisson(mean_numstars) 

195 

196 # In case the Healpix pixel is empty 

197 if nstars == 0: 

198 continue 

199 

200 # Reseed random library 

201 np.random.seed(par.seed + pixel + par.star_dist_seed_offset) 

202 

203 # Positions 

204 x, y = sample_position_uniform(nstars, w, pixel, par.nside_sampling) 

205 ctx.stars.x = np.append(ctx.stars.x, x) 

206 ctx.stars.y = np.append(ctx.stars.y, y) 

207 

208 # Set magnitudes according to catalog 

209 ind = np.random.choice( 

210 np.arange(len(cat_stars_select)), size=nstars, replace=False 

211 ) 

212 

213 # Set magnitudes according to catalog 

214 for f in par.filters: 

215 ctx.stars.magnitude_dict[f] = np.append( 

216 ctx.stars.magnitude_dict[f], cat_stars_select[f][ind] 

217 ) 

218 

219 

220def transform_from_sdss_to_gaia_colours(cat): 

221 # https://www.aanda.org/articles/aa/pdf/2018/08/aa32756-18.pdf Table A.2 

222 colours = np.array( 

223 [ 

224 np.ones(len(cat)), 

225 (cat["g"] - cat["i"]), 

226 (cat["g"] - cat["i"]) ** 2, 

227 (cat["g"] - cat["i"]) ** 3, 

228 ] 

229 ).T 

230 

231 coeffs = np.array([-0.074189, -0.51409, -0.080607, 0.0016001])[np.newaxis, :] 

232 

233 G_gaia_minus_g_sdss = np.sum((colours * coeffs), axis=1) 

234 

235 gaia_G = cat["g"] + G_gaia_minus_g_sdss 

236 

237 return gaia_G 

238 

239 

240def get_star_cat_besancon_gaia_splice(ctx): 

241 # load Besancon map 

242 

243 par = ctx.parameters 

244 w = coordinate_util.wcs_from_parameters(par) 

245 pixels = coordinate_util.get_healpix_pixels( 

246 par.nside_sampling, w, par.size_x, par.size_y 

247 ) 

248 pixarea = hp.nside2pixarea(par.nside_sampling, degrees=True) 

249 

250 # get the scaling for star counts 

251 if hasattr(par, "stars_counts_scale"): 251 ↛ 252line 251 didn't jump to line 252 because the condition on line 251 was never true

252 stars_counts_scale = par.stars_counts_scale 

253 else: 

254 stars_counts_scale = 1 

255 

256 # load besancon 

257 dict_besancon_info = load_besancon_map_info(ctx) 

258 list_ipix_besancon = [] 

259 for pixel in pixels: 

260 theta, phi = hp.pix2ang(par.nside_sampling, pixel) 

261 ipix_besancon = hp.ang2pix(dict_besancon_info["nside"], theta, phi) 

262 list_ipix_besancon.append(ipix_besancon) 

263 

264 besancon_pixels = np.unique(list_ipix_besancon) 

265 dict_cat_besancon = load_besancon_map_pixels(ctx, besancon_pixels) 

266 for key in dict_cat_besancon: 

267 cat = dict_cat_besancon[key] 

268 select = (cat[par.reference_band] > par.stars_mag_min) & ( 

269 cat[par.reference_band] < par.stars_mag_max 

270 ) 

271 cat = cat[select] 

272 

273 cat = at.ensure_cols( 

274 cat, names=["index:i4", "index_gaia:i4", "pos_x", "pos_y", "gaia_G"] 

275 ) 

276 cat["index"] = np.arange(len(cat)) 

277 cat["gaia_G"] = transform_from_sdss_to_gaia_colours(cat) 

278 dict_cat_besancon[key] = cat 

279 

280 # load GAIA 

281 with h5py.File(par.filepath_gaia, "r") as fh5: 

282 cat_gaia = np.array(fh5["data"]) 

283 

284 # get image coordinates for GAIA stars 

285 cat_gaia = at.ensure_cols( 

286 cat_gaia, names=["index:i4", "index_besancon:i4", "x", "y"] 

287 ) 

288 

289 wcs = coordinate_util.tile_in_skycoords( 

290 pixscale=par.pixscale, 

291 ra0=par.ra0, 

292 dec0=par.dec0, 

293 crpix_ra=par.crpix_ra, 

294 crpix_dec=par.crpix_dec, 

295 ) 

296 

297 cat_gaia["x"], cat_gaia["y"] = wcs.all_world2pix(cat_gaia["ra"], cat_gaia["dec"], 0) 

298 cat_gaia["index"] = np.arange(len(cat_gaia)) 

299 

300 # get areas for both catalogues 

301 besancon_area = dict_besancon_info["simulation_area"] 

302 

303 # calculate healpix pixels for Gaia stars 

304 

305 gaia_hppix = coordinate_util.radec2pix( 

306 ra=cat_gaia["ra"], dec=cat_gaia["dec"], nside=par.nside_sampling 

307 ) 

308 

309 list_pixel_besancon = [] 

310 

311 # Loop over pixels 

312 for pixel in pixels: 

313 # get the pixel index in the besancon catalogue 

314 

315 theta, phi = hp.pix2ang(par.nside_sampling, pixel) 

316 ipix_besancon = hp.ang2pix(dict_besancon_info["nside"], theta, phi) 

317 cat_besancon = dict_cat_besancon[f"{ipix_besancon:04d}"] 

318 

319 # current number of stars according to besancon 

320 

321 n_stars_current = np.int32( 

322 np.around(len(cat_besancon) * pixarea / besancon_area) 

323 ) 

324 np.random.seed(par.seed + pixel + par.star_num_seed_offset) 

325 stars_counts_scale = ( 

326 par.stars_counts_scale if hasattr(par, "stars_counts_scale") else 1.0 

327 ) 

328 n_stars_current_noise = np.random.poisson(n_stars_current * stars_counts_scale) 

329 indices_random = np.random.choice( 

330 len(cat_besancon), size=n_stars_current_noise, replace=True 

331 ) 

332 current_besancon = cat_besancon[indices_random].copy() 

333 current_besancon = np.sort(current_besancon, order=["g"]) 

334 np.random.seed(par.seed + pixel + par.star_dist_seed_offset) 

335 x, y = sample_position_uniform( 

336 n_stars_current_noise, w, pixel, par.nside_sampling 

337 ) 

338 current_besancon["index_gaia"] = -99 

339 current_besancon["pos_x"] = x 

340 current_besancon["pos_y"] = y 

341 

342 # select gaia stars in current pixel 

343 

344 current_gaia = cat_gaia[gaia_hppix == pixel].copy() 

345 n_current_gaia = len(current_gaia) 

346 current_besancon_nogaia = current_besancon[n_current_gaia:] 

347 

348 # match gaia and besancon 

349 indices_match_besancon = np.argmin( 

350 ( 

351 ( 

352 current_gaia["phot_g_mean_mag"][:, np.newaxis] 

353 - cat_besancon["gaia_G"][np.newaxis] 

354 ) 

355 ** 2 

356 ), 

357 axis=1, 

358 ) 

359 current_besancon_gaia = cat_besancon[indices_match_besancon].copy() 

360 current_besancon_gaia["index_gaia"] = current_gaia["index"] 

361 current_besancon_gaia["pos_x"] = current_gaia["x"] 

362 current_besancon_gaia["pos_y"] = current_gaia["y"] 

363 

364 current_join_cat = np.concatenate( 

365 [current_besancon_gaia, current_besancon_nogaia] 

366 ) 

367 list_pixel_besancon.append(current_join_cat) 

368 

369 cat_full = np.concatenate(list_pixel_besancon) 

370 ctx.stars.x = np.append(ctx.stars.x, cat_full["pos_x"]) 

371 ctx.stars.y = np.append(ctx.stars.y, cat_full["pos_y"]) 

372 

373 # Set magnitudes according to catalog 

374 for f in par.filters: 

375 ctx.stars.magnitude_dict[f] = np.append( 

376 ctx.stars.magnitude_dict[f], cat_full[f] 

377 ) 

378 

379 

380STAR_GENERATOR = { 

381 "besancon_map": get_star_cat_besancon, 

382 "besancon_gaia_splice": get_star_cat_besancon_gaia_splice, 

383} 

384 

385 

386class Plugin(BasePlugin): 

387 """ 

388 Draw stars for the r-band from a catalog created using the Besancon model of the 

389 Milky Way. 

390 

391 :param besancon_cat_name: Path to the catalog of stars drawn from the model stored 

392 as fits-file. The catalog must contain a column called 'r' giving the magnitude 

393 of the stars in the r-band. Furthermore, the header of the HDU in which the 

394 catalog is stored must contain a field labeled 'AREA' which states the area (in 

395 sq. deg) the catalog covers. 

396 :param seed: General seed. 

397 :param star_num_seed_offset: Seed offset before number of stars is drawn. 

398 :param stars_mag_min: Minimum magnitude of stars in the r-band. 

399 :param stars_mag_max: Maximum magnitude of stars in the r-band. 

400 :param star_dist_seed_offset: Seed offset before positions of stars are drawn. 

401 :param star_nphot_seed_offset: Seed offset before numbers of photons are calculated 

402 for stars. 

403 :param n_exp: Number of single exposures in coadded image. 

404 :param magzero: Magnitude zeropoint. 

405 :param gain: Gain in electrons per ADU. 

406 :param exp_time_file_name: File name of an exposure time map. 

407 :param size_x: Size of the image in x-direction. 

408 :param size_y: Size of the image in y-direction. 

409 :param maps_remote_dir: Remote directory where maps are stored. 

410 """ 

411 

412 def __call__(self): 

413 par = self.ctx.parameters 

414 

415 if not hasattr(par, "crpix_ra"): 

416 par.crpix_ra = par.size_x / 2 + 0.5 

417 if not hasattr(par, "crpix_dec"): 

418 par.crpix_dec = par.size_y / 2 + 0.5 

419 

420 # Initialize star catalog 

421 self.ctx.stars = sampling_util.Catalog() 

422 self.ctx.stars.columns = ["id", "x", "y"] 

423 

424 self.ctx.stars.x = np.zeros(0, dtype=np.float32) 

425 self.ctx.stars.y = np.zeros(0, dtype=np.float32) 

426 self.ctx.stars.magnitude_dict = {} 

427 for f in par.filters: 

428 self.ctx.stars.magnitude_dict[f] = np.zeros(0, dtype=np.float32) 

429 

430 # Get magnitudes and position 

431 generator = STAR_GENERATOR[par.star_catalogue_type] 

432 generator(self.ctx) 

433 

434 # Mask stars not on image 

435 mask = ( 

436 (self.ctx.stars.x >= 0) 

437 & (self.ctx.stars.x < par.size_x) 

438 & (self.ctx.stars.y >= 0) 

439 & (self.ctx.stars.y < par.size_y) 

440 ) 

441 

442 self.ctx.numstars = np.sum(mask) 

443 self.ctx.stars.id = np.arange(self.ctx.numstars) 

444 

445 self.ctx.stars.x = self.ctx.stars.x[mask] 

446 self.ctx.stars.y = self.ctx.stars.y[mask] 

447 

448 for f in par.filters: 

449 self.ctx.stars.magnitude_dict[f] = self.ctx.stars.magnitude_dict[f][mask] 

450 

451 def __str__(self): 

452 return NAME