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-17 08:48 +0000
1# Copyright (C) 2016 ETH Zurich, Institute for Astronomy
3"""
4Created on Jul 20, 2016
5author: Tomasz Kacprzak
6"""
8import warnings
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
17from .. import coordinate_util, io_util, sampling_util
18from ..sampling_util import sample_position_uniform
20if six.PY2: # pragma: no cover
21 import cPickle as pickle
22else:
23 import pickle
26NAME = "draw stars"
29def get_interp_nearest(ra_new, dec_new, dict_besancon_info):
30 """
32 :param ra_new:
33 :param dec_new:
34 :param dict_besancon_info:
35 :return:
36 """
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
49 closest_cat = dict_besancon_info["healpix_list"][ipix]
51 return closest_cat
54def load_besancon_map_info(ctx):
55 par = ctx.parameters
57 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir)
59 dict_besancon_info = {}
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"]))
65 return dict_besancon_info
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:
73 dict_besancon_info={'simulation_area': area_store, 'nside': nside,
74 'healpix_list': list_df, 'healpix_mask': hp_map}
76 simulation_area - area in sq. deg corresponding to the simulation catalog covers
78 nside - resolution of simulation sampling
80 healpix_list - a list of numpy record arrays, each element corresponds to a pixel on
81 HEALPix grid with nside=nside
83 hp_map - HEALPix map showing the coverage of simulations, ring=False
84 """
86 par = ctx.parameters
88 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir)
90 dist_besancon_pix_cat = {}
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}"])
99 dist_besancon_pix_cat[f"{ipix:04d}"] = besancon_pix_cat
101 return dist_besancon_pix_cat
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:
109 dict_besancon_info = {'simulation_area': area_store, 'nside': nside,
110 'healpix_list': list_df, 'healpix_mask': hp_map}
112 simulation_area - area in sq. deg corresponding to the simulation catalog covers
114 nside - resolution of simulation sampling
116 healpix_list - a list of numpy record arrays, each element corresponds
117 to a pixel on HEALPix grid with nside=nside
119 hp_map - HEALPix map showing the coverage of simulations, ring=False
120 """
122 par = ctx.parameters
124 filepath_map = io_util.get_abs_path(par.besancon_map_path, par.maps_remote_dir)
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)
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 )
148 return dict_besancon_info
151def get_star_cat_besancon(ctx):
152 par = ctx.parameters
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 )
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)
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"]
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 )
179 cat_stars_select = cat_stars[mask]
181 max_numstars = np.count_nonzero(mask)
183 # Mean number of stars within a Healpix pixel
184 mean_numstars = np.int32(np.around(max_numstars * pixarea / simulation_area))
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)
196 # In case the Healpix pixel is empty
197 if nstars == 0:
198 continue
200 # Reseed random library
201 np.random.seed(par.seed + pixel + par.star_dist_seed_offset)
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)
208 # Set magnitudes according to catalog
209 ind = np.random.choice(
210 np.arange(len(cat_stars_select)), size=nstars, replace=False
211 )
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 )
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
231 coeffs = np.array([-0.074189, -0.51409, -0.080607, 0.0016001])[np.newaxis, :]
233 G_gaia_minus_g_sdss = np.sum((colours * coeffs), axis=1)
235 gaia_G = cat["g"] + G_gaia_minus_g_sdss
237 return gaia_G
240def get_star_cat_besancon_gaia_splice(ctx):
241 # load Besancon map
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)
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
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)
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]
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
280 # load GAIA
281 with h5py.File(par.filepath_gaia, "r") as fh5:
282 cat_gaia = np.array(fh5["data"])
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 )
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 )
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))
300 # get areas for both catalogues
301 besancon_area = dict_besancon_info["simulation_area"]
303 # calculate healpix pixels for Gaia stars
305 gaia_hppix = coordinate_util.radec2pix(
306 ra=cat_gaia["ra"], dec=cat_gaia["dec"], nside=par.nside_sampling
307 )
309 list_pixel_besancon = []
311 # Loop over pixels
312 for pixel in pixels:
313 # get the pixel index in the besancon catalogue
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}"]
319 # current number of stars according to besancon
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
342 # select gaia stars in current pixel
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:]
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"]
364 current_join_cat = np.concatenate(
365 [current_besancon_gaia, current_besancon_nogaia]
366 )
367 list_pixel_besancon.append(current_join_cat)
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"])
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 )
380STAR_GENERATOR = {
381 "besancon_map": get_star_cat_besancon,
382 "besancon_gaia_splice": get_star_cat_besancon_gaia_splice,
383}
386class Plugin(BasePlugin):
387 """
388 Draw stars for the r-band from a catalog created using the Besancon model of the
389 Milky Way.
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 """
412 def __call__(self):
413 par = self.ctx.parameters
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
420 # Initialize star catalog
421 self.ctx.stars = sampling_util.Catalog()
422 self.ctx.stars.columns = ["id", "x", "y"]
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)
430 # Get magnitudes and position
431 generator = STAR_GENERATOR[par.star_catalogue_type]
432 generator(self.ctx)
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 )
442 self.ctx.numstars = np.sum(mask)
443 self.ctx.stars.id = np.arange(self.ctx.numstars)
445 self.ctx.stars.x = self.ctx.stars.x[mask]
446 self.ctx.stars.y = self.ctx.stars.y[mask]
448 for f in par.filters:
449 self.ctx.stars.magnitude_dict[f] = self.ctx.stars.magnitude_dict[f][mask]
451 def __str__(self):
452 return NAME