Coverage for src/ufig/plugins/add_psf.py: 96%
294 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-12 19:08 +0000
1# Copyright (c) 2015 ETH Zurich, Institute of Astronomy, Claudio Bruderer
2# <claudio.bruderer@phys.ethz.ch>
4"""
5Created on Mar 20, 2015
6@author: Claudio Bruderer
7adapted by Silvan Fischbacher, 2024
8"""
10import h5py
11import healpy as hp
12import numba as nb
13import numpy as np
14from cosmic_toolbox import logger
15from ivy.plugin.base_plugin import BasePlugin
16from scipy import optimize
18from ufig import io_util
19from ufig.psf_estimation import correct_brighter_fatter
21from .. import coordinate_util
23LOGGER = logger.get_logger(__file__)
25NAME = "add psf"
28@nb.njit(nb.types.int32[:](nb.types.float32[:, :], nb.types.float32[:, :]), cache=True)
29def numba_min_dist(X, Y):
30 """
31 Finds the index of the closest point in Y for each point in X
32 """
33 m = X.shape[0]
34 D = np.zeros((m,), dtype=np.int32)
35 for i in nb.prange(m):
36 arg_min = np.argmin(np.sum((X[i] - Y) ** 2, axis=1))
37 D[i] = arg_min
39 return D
42def load_psf_skymaps(psf_maps, par):
43 """
44 Load maps of the PSF across the sky.
46 :param psf_maps: File name of the psf maps (0: r50-map, 1: e1-map, 2: e2-map)
47 :param par: ctx.parameters containing potential fudge factors for PSF maps
49 :return r50_map: r50-map containing flux radius (50%)-variations across the survey
50 area
51 :return e1_map: Ellipticity 1-component-map containing variations across the survey
52 area
53 :return e2_map: Ellipticity 2-component-map containing variations across the survey
54 area
55 """
57 maps = io_util.load_hpmap(psf_maps, par.maps_remote_dir)
58 r50_map = par.psf_r50_factor * maps[0].astype(np.float32) + par.psf_r50_shift
59 r50_map[r50_map < 0] = 0.0
60 e1_map = par.psf_e1_prefactor * (
61 par.psf_e1_factor * maps[1].astype(np.float32) + par.psf_e1_shift
62 )
63 e2_map = par.psf_e2_factor * maps[2].astype(np.float32) + par.psf_e2_shift
64 return r50_map, e1_map, e2_map
67def psf_from_sky_maps(r50_map, e1_map, e2_map, psf_beta, w, x, y, psf_flux_ratio):
68 """
69 Evaluate PSF maps at input positions.
71 :param r50_map: r50-map containing flux radius (50%)-variations across the survey
72 area
73 :param e1_map: Ellipticity 1-component-map containing variations across the survey
74 area
75 :param e2_map: Ellipticity 2-component-map containing variations across the survey
76 area
77 :param psf_beta: PSF beta parameter evaluated at every object's position
78 :param w: wcs-object containing all the relevant wcs-transformation information
79 :param x: pixel x-coordinate
80 :param y: pixel y-coordinate
81 :param psf_flux_ratio: Flux ratio of the different PSF Moffat component(s) relative
82 to the total flux
84 :return r50: PSF r50 (flux radius 50%) evaluated at the input positions
85 :return beta: PSF beta parameter evaluated at the input positions
86 :return e1: PSF ellipticity 1-component evaluated at the input positions
87 :return e2: PSF ellipticity 2-component evaluated at the input positions
88 """
90 # +0.5 is to convert it into origin-1 convention
91 theta, phi = coordinate_util.xy2thetaphi(w, x, y)
93 beta = np.full((x.size, len(psf_beta)), psf_beta, np.float32)
95 # PSF interpolation
96 nside = hp.get_nside(r50_map)
98 pix_indices = hp.ang2pix(nside=nside, theta=theta, phi=phi, nest=False)
100 r50 = r50_map[pix_indices]
101 e1 = e1_map[pix_indices]
102 e2 = e2_map[pix_indices]
104 # Convert psf_r50 to psf_fwhm (multiple betas)
105 psf_fwhm = moffat_r502fwhm(r50, psf_beta, psf_flux_ratio)
107 psf_r50_indiv = np.empty((len(psf_beta), x.size), dtype=np.float32)
108 for i in range(len(psf_beta)):
109 psf_r50_indiv[i] = moffat_fwhm2r50(psf_fwhm, psf_beta[i], psf_flux_ratio)
111 return r50, beta, e1, e2, psf_r50_indiv, psf_fwhm
114def get_moffat_maps_psf(obj_name, ctx):
115 """
116 Set the PSF for objects (stars or galaxies) using maps of the PSF across the sky.
118 :param obj_name: Name of the objects the PSF is evaluated for
119 :param ctx: Ivy-context containing the catalog of the object properties
121 :return r50: PSF r50 (flux radius 50%) evaluated at the positions of the input
122 objects
123 :return beta: PSF beta parameter evaluated at the positions of the input objects
124 :return e1: PSF ellipticity 1-component evaluated at the positions of the input
125 objects
126 :return e2: PSF ellipticity 2-component evaluated at the positions of the input
127 objects
128 """
130 par = ctx.parameters
132 obj = getattr(ctx, obj_name)
134 if not isinstance(par.psf_beta, list):
135 par.psf_beta = [par.psf_beta]
137 if len(par.psf_beta) == 1: 137 ↛ 140line 137 didn't jump to line 140 because the condition on line 137 was always true
138 par.psf_flux_ratio = 1
140 if not hasattr(ctx, "psf_r50_map"):
141 ctx.psf_r50_map, ctx.psf_e1_map, ctx.psf_e2_map = load_psf_skymaps(
142 par.psf_maps, par
143 )
145 w = coordinate_util.tile_in_skycoords(
146 pixscale=par.pixscale,
147 ra0=par.ra0,
148 dec0=par.dec0,
149 crpix_ra=par.crpix_ra,
150 crpix_dec=par.crpix_dec,
151 )
153 x = obj.x
154 y = obj.y
156 psf_r50, psf_beta, psf_e1, psf_e2, psf_r50_indiv, psf_fwhm = psf_from_sky_maps(
157 r50_map=ctx.psf_r50_map,
158 e1_map=ctx.psf_e1_map,
159 e2_map=ctx.psf_e2_map,
160 psf_beta=par.psf_beta,
161 w=w,
162 x=x,
163 y=y,
164 psf_flux_ratio=par.psf_flux_ratio,
165 )
166 ctx.psf_column_names = [
167 "psf_r50",
168 "psf_beta",
169 "psf_e1",
170 "psf_e2",
171 "psf_r50_indiv",
172 "psf_fwhm",
173 ]
175 obj.psf_r50 = psf_r50
176 obj.psf_beta = psf_beta
177 obj.psf_e1 = psf_e1
178 obj.psf_e2 = psf_e2
179 obj.psf_r50_indiv = psf_r50_indiv
180 obj.psf_fwhm = psf_fwhm
183def apply_psf_parameter_scalings(psf_par, par):
184 for name in psf_par.dtype.names:
185 if name in par.psf_cnn_factors: 185 ↛ 186line 185 didn't jump to line 186
186 psf_par[name] = (
187 par.psf_cnn_factors[name][0]
188 + par.psf_cnn_factors[name][1] * psf_par[name]
189 )
192def ensure_valid_psf_beta(psf_beta):
193 psf_beta[:] = np.clip(psf_beta, a_min=1.2, a_max=np.inf)
196def ensure_valid_psf_flux_ratio(psf_flux_ratio):
197 psf_flux_ratio[:] = np.clip(psf_flux_ratio, a_min=0.0, a_max=1.0)
198 psf_flux_ratio[:] = np.nan_to_num(psf_flux_ratio, nan=1.0)
201def ensure_valid_psf_fwhm(psf_fwhm):
202 select = psf_fwhm <= 0
203 psf_fwhm[select] = 0.0001
206def ensure_valid_psf_ellip(psf_e1, psf_e2):
207 psf_ellip = np.sqrt(psf_e1**2 + psf_e2**2)
208 select = psf_ellip > 1
209 psf_e1[select] /= 1.01 * psf_ellip[select]
210 psf_e2[select] /= 1.01 * psf_ellip[select]
213def get_moffat_coadd_psf_cnn(obj_name, ctx):
214 # Parameters
215 par = ctx.parameters
217 # Objects
218 obj = getattr(ctx, obj_name)
220 # Predict PSF at object position
221 predict_args = dict(
222 position_xy=np.stack((obj.x, obj.y), axis=-1),
223 filepath_psfmodel=io_util.get_abs_path(
224 par.filepath_psfmodel_input, par.maps_remote_dir
225 ),
226 )
228 if par.psf_type == "coadd_moffat_cnn_robust" or par.psf_type == "coadd_moffat_cnn":
229 if par.psf_type == "coadd_moffat_cnn_robust": # pragma: no cover
230 LOGGER.warning("Using deprecated PSF type name 'coadd_moffat_cnn_robust'")
231 LOGGER.warning("Please use 'coadd_moffat_cnn' instead")
233 from ufig.psf_estimation import psf_estimation_coadd_cnn as psf_est
235 col_ending = "_ipt"
237 psf_par_pred = next(psf_est.predict_psf_with_file(**predict_args))[0]
239 # Scale predictions
240 apply_psf_parameter_scalings(psf_par_pred, par)
242 # Set PSF beta
243 n_beta = len(
244 list(
245 filter(
246 lambda col: col.startswith("psf_beta_") and col.endswith(col_ending),
247 psf_par_pred.dtype.names,
248 )
249 )
250 )
252 if n_beta > 0: 252 ↛ 253line 252 didn't jump to line 253 because the condition on line 252 was never true
253 obj.psf_beta = np.empty((len(psf_par_pred), n_beta), dtype=np.float32)
255 for i_beta in range(n_beta):
256 obj.psf_beta[:, i_beta] = psf_par_pred[f"psf_beta_{i_beta + 1}{col_ending}"]
258 else:
259 obj.psf_beta = np.full(
260 (len(psf_par_pred), len(par.psf_beta)), par.psf_beta, dtype=np.float32
261 )
263 # Set other PSF parameters
264 psf_par_names = [
265 "psf_flux_ratio",
266 "psf_fwhm",
267 "psf_e1",
268 "psf_e2",
269 "psf_f1",
270 "psf_f2",
271 "psf_g1",
272 "psf_g2",
273 "psf_kurtosis",
274 ]
276 for psf_par_name in psf_par_names:
277 psf_par_name_ipt = psf_par_name + col_ending
279 if psf_par_name_ipt in psf_par_pred.dtype.names:
280 psf_par_col = psf_par_pred[psf_par_name_ipt]
282 else:
283 psf_par_col = np.full(
284 len(psf_par_pred),
285 getattr(ctx.parameters, psf_par_name),
286 dtype=np.float32,
287 )
289 setattr(obj, psf_par_name, psf_par_col)
291 # Ensure valid beta, flux ratio, size and ellipticity
292 ensure_valid_psf_beta(obj.psf_beta)
293 ensure_valid_psf_flux_ratio(obj.psf_flux_ratio)
294 ensure_valid_psf_fwhm(obj.psf_fwhm)
295 ensure_valid_psf_ellip(obj.psf_e1, obj.psf_e2)
297 # Set half-light radii and offsets between profiles
298 obj.psf_r50_indiv = (
299 obj.psf_fwhm[:, np.newaxis]
300 * np.sqrt((2 ** (1 / (obj.psf_beta - 1)) - 1) / (2 ** (1 / obj.psf_beta) - 1))
301 / 2
302 )
303 obj.psf_r50_indiv = obj.psf_r50_indiv.T
304 obj.psf_dx_offset = np.zeros_like(obj.psf_r50_indiv).T
305 obj.psf_dy_offset = np.zeros_like(obj.psf_r50_indiv).T
307 # Apply brighter-fatter corrections
308 if ctx.parameters.psfmodel_corr_brighter_fatter is not None: 308 ↛ 334line 308 didn't jump to line 334
309 apply_brighter_fatter = False
311 if obj_name == "stars":
312 apply_brighter_fatter = True
314 if (obj_name == "galaxies") and (
315 "apply_to_galaxies" in ctx.parameters.psfmodel_corr_brighter_fatter
316 ):
317 apply_brighter_fatter = ctx.parameters.psfmodel_corr_brighter_fatter[
318 "apply_to_galaxies"
319 ]
321 if apply_brighter_fatter:
322 (
323 obj.psf_fwhm,
324 obj.psf_e1,
325 obj.psf_e2,
326 ) = correct_brighter_fatter.brighter_fatter_add(
327 col_mag=obj.mag,
328 col_fwhm=obj.psf_fwhm,
329 col_e1=obj.psf_e1,
330 col_e2=obj.psf_e2,
331 dict_corr=ctx.parameters.psfmodel_corr_brighter_fatter,
332 )
334 ctx.psf_column_names = [
335 "psf_flux_ratio",
336 "psf_fwhm",
337 "psf_e1",
338 "psf_e2",
339 "psf_f1",
340 "psf_f2",
341 "psf_g1",
342 "psf_g2",
343 "psf_kurtosis",
344 "psf_r50_indiv",
345 "psf_dy_offset",
346 "psf_dx_offset",
347 ]
350def get_moffat_coadd_psf_cnn_from_file(obj_name, ctx):
351 # Parameters
352 par = ctx.parameters
354 # Objects
355 obj = getattr(ctx, obj_name)
357 # Predict PSF at object position
358 predict_args = dict(
359 position_xy=np.stack((obj.x, obj.y), axis=-1),
360 filepath_psfmodel=par.filepath_psfmodel_input,
361 )
363 col_ending = "_ipt"
365 psfcat = h5py.File(predict_args["filepath_psfmodel"])
367 xy_ucat = np.stack((obj.x, obj.y)).T
368 xy_ucat += 0.5
369 xy_psf = np.stack((psfcat["grid_psf"]["X_IMAGE"], psfcat["grid_psf"]["Y_IMAGE"])).T
370 min_dist = numba_min_dist(xy_ucat.astype(np.float32), xy_psf.astype(np.float32))
371 psf_params = [
372 "psf_flux_ratio_ipt",
373 "psf_fwhm_ipt",
374 "psf_e1_ipt",
375 "psf_e2_ipt",
376 "psf_f1_ipt",
377 "psf_f2_ipt",
378 "psf_g1_ipt",
379 "psf_g2_ipt",
380 ]
381 psf_par_pred = psfcat["grid_psf"][:][psf_params][min_dist]
383 # Set PSF beta
384 n_beta = len(
385 list(
386 filter(
387 lambda col: col.startswith("psf_beta_") and col.endswith(col_ending),
388 psf_par_pred.dtype.names,
389 )
390 )
391 )
393 if n_beta > 0: 393 ↛ 394line 393 didn't jump to line 394 because the condition on line 393 was never true
394 obj.psf_beta = np.empty((len(psf_par_pred), n_beta), dtype=np.float32)
396 for i_beta in range(n_beta):
397 obj.psf_beta[:, i_beta] = psf_par_pred[f"psf_beta_{i_beta + 1}{col_ending}"]
399 else:
400 obj.psf_beta = np.full(
401 (len(psf_par_pred), len(par.psf_beta)), par.psf_beta, dtype=np.float32
402 )
404 # Set other PSF parameters
405 psf_par_names = [
406 "psf_flux_ratio",
407 "psf_fwhm",
408 "psf_e1",
409 "psf_e2",
410 "psf_f1",
411 "psf_f2",
412 "psf_g1",
413 "psf_g2",
414 "psf_kurtosis",
415 ]
417 for psf_par_name in psf_par_names:
418 psf_par_name_ipt = psf_par_name + col_ending
420 if psf_par_name_ipt in psf_par_pred.dtype.names:
421 psf_par_col = psf_par_pred[psf_par_name_ipt]
423 else:
424 psf_par_col = np.full(
425 len(psf_par_pred),
426 getattr(ctx.parameters, psf_par_name),
427 dtype=np.float32,
428 )
430 setattr(obj, psf_par_name, psf_par_col)
432 # Ensure valid beta, flux ratio, size and ellipticity
433 ensure_valid_psf_beta(obj.psf_beta)
434 ensure_valid_psf_flux_ratio(obj.psf_flux_ratio)
435 ensure_valid_psf_fwhm(obj.psf_fwhm)
436 ensure_valid_psf_ellip(obj.psf_e1, obj.psf_e2)
438 # Set half-light radii and offsets between profiles
439 obj.psf_r50_indiv = (
440 obj.psf_fwhm[:, np.newaxis]
441 * np.sqrt((2 ** (1 / (obj.psf_beta - 1)) - 1) / (2 ** (1 / obj.psf_beta) - 1))
442 / 2
443 )
444 obj.psf_r50_indiv = obj.psf_r50_indiv.T
445 obj.psf_dx_offset = np.zeros_like(obj.psf_r50_indiv).T
446 obj.psf_dy_offset = np.zeros_like(obj.psf_r50_indiv).T
448 # Apply brighter-fatter corrections
449 if ctx.parameters.psfmodel_corr_brighter_fatter is not None: 449 ↛ 475line 449 didn't jump to line 475
450 apply_brighter_fatter = False
452 if obj_name == "stars":
453 apply_brighter_fatter = True
455 if (obj_name == "galaxies") and (
456 "apply_to_galaxies" in ctx.parameters.psfmodel_corr_brighter_fatter
457 ):
458 apply_brighter_fatter = ctx.parameters.psfmodel_corr_brighter_fatter[
459 "apply_to_galaxies"
460 ]
462 if apply_brighter_fatter:
463 (
464 obj.psf_fwhm,
465 obj.psf_e1,
466 obj.psf_e2,
467 ) = correct_brighter_fatter.brighter_fatter_add(
468 col_mag=obj.mag,
469 col_fwhm=obj.psf_fwhm,
470 col_e1=obj.psf_e1,
471 col_e2=obj.psf_e2,
472 dict_corr=ctx.parameters.psfmodel_corr_brighter_fatter,
473 )
475 ctx.psf_column_names = [
476 "psf_flux_ratio",
477 "psf_fwhm",
478 "psf_e1",
479 "psf_e2",
480 "psf_f1",
481 "psf_f2",
482 "psf_g1",
483 "psf_g2",
484 "psf_kurtosis",
485 "psf_r50_indiv",
486 "psf_dy_offset",
487 "psf_dx_offset",
488 ]
491def sample_psf_moffat_constant(obj_name, ctx):
492 """
493 Evaluate a constant, Moffat-PSF field at the positions of different objects (stars
494 or galaxies).
496 :param obj: Name of the objects the PSF is evaluated for
497 :param ctx: Ivy-context containing the catalog of the object properties
498 :return r50: PSF r50 (flux radius 50%) evaluated at the positions of the input
499 objects
500 :return beta: PSF beta parameter evaluated at the positions of the input objects
501 :return e1: PSF ellipticity 1-component evaluated at the positions of the input
502 objects
503 :return e2: PSF ellipticity 2-component evaluated at the positions of the input
504 objects
505 """
507 par = ctx.parameters
509 # Objects
510 obj = getattr(ctx, obj_name)
512 if not isinstance(par.psf_beta, list):
513 par.psf_beta = [par.psf_beta]
515 if len(par.psf_beta) == 1:
516 par.psf_flux_ratio = 1
518 numobj = obj.x.size
520 psf_fwhm = np.full(numobj, par.seeing / par.pixscale, dtype=np.float32)
521 psf_r50 = moffat_fwhm2r50(psf_fwhm, par.psf_beta, par.psf_flux_ratio)
522 psf_beta = np.full((numobj, len(par.psf_beta)), par.psf_beta, dtype=np.float32)
523 psf_e1 = np.full(numobj, par.psf_e1, dtype=np.float32)
524 psf_e2 = np.full(numobj, par.psf_e2, dtype=np.float32)
526 psf_r50_indiv = np.empty((len(par.psf_beta), numobj))
527 for i in range(len(par.psf_beta)):
528 psf_r50_indiv[i] = moffat_fwhm2r50(
529 psf_fwhm, par.psf_beta[i], par.psf_flux_ratio
530 )
532 obj.psf_dx_offset = np.zeros([len(obj.x), 2])
533 obj.psf_dy_offset = np.zeros([len(obj.y), 2])
535 ctx.psf_column_names = [
536 "psf_r50",
537 "psf_beta",
538 "psf_e1",
539 "psf_e2",
540 "psf_r50_indiv",
541 "psf_fwhm",
542 ]
544 obj.psf_r50 = psf_r50
545 obj.psf_beta = psf_beta
546 obj.psf_e1 = psf_e1
547 obj.psf_e2 = psf_e2
548 obj.psf_r50_indiv = psf_r50_indiv
549 obj.psf_fwhm = psf_fwhm
552def moffat_r502fwhm(psf_r50, psf_beta, psf_flux_ratio):
553 """
554 Computes the FWHM for a given r50 and beta parameter for the PSF assuming a Moffat
555 distribution
557 :param psf_r50: r50 of the PSF in units of pixels
558 :param psf_beta: beta parameter for a Moffat distribution
559 :return: FWHM of the PSF (assuming a Moffat distribution) in units of pixels
560 """
562 def convert(r50, beta):
563 return 2 * r50 * np.sqrt((2 ** (1 / beta) - 1) / (2 ** (1 / (beta - 1)) - 1))
565 if isinstance(psf_beta, np.ndarray):
566 psf_fwhm = convert(psf_r50, psf_beta)
568 else:
569 if not isinstance(psf_beta, list):
570 psf_beta = [psf_beta]
572 if len(psf_beta) == 1:
573 psf_fwhm = convert(psf_r50, psf_beta[0])
574 else:
575 psf_fwhm = multiple_moffat_fwhm(psf_r50, psf_beta, psf_flux_ratio)
577 return psf_fwhm
580def moffat_fwhm2r50(psf_fwhm, psf_beta, psf_flux_ratio):
581 """
582 Computes the r50 for a given FWHM and beta parameter for the PSF assuming a Moffat
583 distribution
585 :param psf_fwhm: FWHM of the PSF in units of pixels
586 :param psf_beta: beta parameter for a Moffat distribution
587 :return: r50 of the PSF (assuming a Moffat distribution) in units of pixels
588 """
590 def convert(fwhm, beta):
591 return fwhm / 2 * np.sqrt((2 ** (1 / (beta - 1)) - 1) / (2 ** (1 / beta) - 1))
593 if isinstance(psf_beta, np.ndarray):
594 psf_r50 = convert(psf_fwhm, psf_beta)
596 else:
597 if not isinstance(psf_beta, list):
598 psf_beta = [psf_beta]
600 if len(psf_beta) == 1:
601 psf_r50 = convert(psf_fwhm, psf_beta[0])
602 else:
603 psf_r50 = multiple_moffat_r50(psf_fwhm, psf_beta, psf_flux_ratio)
605 return psf_r50
608def moffat_fwhm2alpha(psf_fwhm, psf_beta):
609 """
610 Computes the alpha parameter for a given FWHM and beta parameter for a Moffat
611 distribution
613 :param psf_fwhm: FWHM of the PSF in units of pixels
614 :param psf_beta: beta parameter for a Moffat distribution
615 :return: alpha of the Moffat profile in units of pixels
616 """
618 if isinstance(psf_beta, list):
619 psf_beta = psf_beta[0]
621 alpha = psf_fwhm / 2 / np.sqrt(2 ** (1 / psf_beta) - 1)
622 return alpha
625def moffat_r502alpha(psf_r50, psf_beta):
626 """
627 Computes the alpha parameter for a given r50 and beta parameter for a Moffat
628 distribution
630 :param psf_r50: r50 of the PSF in units of pixels
631 :param psf_beta: beta parameter for a Moffat distribution
632 :return: alpha of the Moffat profile in units of pixels
633 """
635 if isinstance(psf_beta, list):
636 psf_beta = psf_beta[0]
638 alpha = psf_r50 / np.sqrt(2 ** (1 / (psf_beta - 1)) - 1)
639 return alpha
642def moffat_alpha2fwhm(psf_alpha, psf_beta):
643 """
644 Computes the fwhm (in the units alpha is in) for a given alpha and beta parameter
645 for a Moffat distribution
647 :param psf_alpha: alpha of the PSF in units of pixels
648 :param psf_beta: beta parameter for a Moffat distribution
649 :return: alpha of the Moffat profile in units of pixels
650 """
652 if isinstance(psf_beta, list):
653 psf_beta = psf_beta[0]
655 psf_fwhm = 2 * psf_alpha * np.sqrt(2 ** (1 / psf_beta) - 1)
657 return psf_fwhm
660def moffat_alpha2r50(psf_alpha, psf_beta):
661 """
662 Computes the r50 (in units of alpha) for a given alpha and beta parameter for a
663 Moffat distribution
665 :param psf_alpha: alpha of the PSF in units of pixels
666 :param psf_beta: beta parameter for a Moffat distribution
667 :return: alpha of the Moffat profile in units of pixels
668 """
670 if isinstance(psf_beta, list):
671 psf_beta = psf_beta[0]
673 psf_r50 = psf_alpha * np.sqrt(2 ** (1 / (psf_beta - 1)) - 1)
675 return psf_r50
678def moffat_profile_integrated(r, psf_beta, psf_flux_ratio):
679 """
680 Evaluates the integrated "1 - Moffat profile" within a radius r (in units of Moffat
681 alpha), which potentially can consist of multiple components, and returns the flux
682 in units of the total flux.
684 :param x: Radius in units of Moffat alpha
685 :param psf_beta: Moffat beta parameter(s)
686 :param psf_flux_ratio: Fraction of flux in the first component
687 :return: Flux within a radius r in units of the total flux
688 """
690 if not isinstance(psf_beta, list):
691 psf_beta = [psf_beta]
693 if len(psf_beta) == 1:
694 psf_flux_ratio = 1
696 value = 0.0
697 for i in range(len(psf_beta)):
698 value += psf_flux_ratio / (1 + r**2) ** (psf_beta[i] - 1)
699 psf_flux_ratio = 1 - psf_flux_ratio
701 return value
704def _moffat_subtracted(x, ppsf_beta, ppsf_flux_ratio):
705 return moffat_profile_integrated(x, ppsf_beta, ppsf_flux_ratio) - 0.5
708def multiple_moffat_r50(psf_fwhm, psf_beta, psf_flux_ratio):
709 """
710 Solve the nonlinear equation to recover an effective r50 of a PSF profile
711 consisting of potentially multiple Moffat profiles.
713 :param psf_fwhm: FWHM of the Moffat profile in pixels
714 :param psf_beta: Moffat beta parameter(s)
715 :param psf_flux_ratio: Fraction of flux in the first component
716 :return: r50 of the total PSF profile
717 """
719 r50 = optimize.brentq(
720 _moffat_subtracted, 0.0, 10 * np.mean(psf_fwhm), args=(psf_beta, psf_flux_ratio)
721 )
723 alpha = moffat_fwhm2alpha(psf_fwhm, psf_beta[0])
724 r50 = r50 * alpha
726 return r50
729def multiple_moffat_fwhm(psf_r50, psf_beta, psf_flux_ratio):
730 """
731 Solve the nonlinear equation to recover an effective fwhm of a PSF profile
732 consisting of potentially multiple Moffat profiles.
733 This function is similar to ufig.plugins.add_psf.multiple_moffat_r50().
735 :param psf_fwhm: FWHM of the Moffat profile in pixels
736 :param psf_beta: Moffat beta parameter(s)
737 :param psf_flux_ratio: Fraction of flux in the first component
738 :return: r50 of the total PSF profile
739 """
741 r50_reduced = optimize.brentq(
742 _moffat_subtracted, 0.0, 10 * np.mean(psf_r50), args=(psf_beta, psf_flux_ratio)
743 )
745 alpha = psf_r50 / r50_reduced
746 fwhm = moffat_alpha2fwhm(alpha, psf_beta[0])
748 return fwhm
751def update_psf_for_current_filter(obj, ctx):
752 for col in ctx.psf_column_names:
753 col_bands = f"{col}_dict"
754 if not hasattr(obj, col_bands):
755 col_dict = setattr(obj, col_bands, {})
757 col_dict = getattr(obj, col_bands)
758 col_dict[ctx.current_filter] = getattr(obj, col)
761PSF_GENERATOR = {
762 "constant_moffat": sample_psf_moffat_constant,
763 "maps_moffat": get_moffat_maps_psf,
764 "coadd_moffat_cnn": get_moffat_coadd_psf_cnn,
765 "coadd_moffat_cnn_robust": get_moffat_coadd_psf_cnn, # backwards compatibility
766 "coadd_moffat_cnn_read": get_moffat_coadd_psf_cnn_from_file, # keep
767}
770class Plugin(BasePlugin):
771 """
772 The PSF is applied to the input galaxy and star catalogs by evaluating the PSF size
773 and ellipticity fields at every objects position
775 :param psf_type: whether a constant or variable psf is added
776 :param psf_beta: beta parameter for Moffat PSF profile
777 :param seeing: seeing (arcsec)
778 :param psf_e1: e1 for the PSF
779 :param psf_e2: e2 for the PSF
780 :param psf_maps: File name of the psf maps (0: r50-map, 1: e1-map, 2: e2-map)
781 :param maps_remote_dir: Remote directory used in case the file name of the psf maps
782 is not absolute
783 :param pixscale: pixel scale (arcsec/pixel)
784 :param size_x: size of image in x-direction (pixel)
785 :param size_y: size of image in y-direction (pixel)
786 :param ra0: right ascension at the center of the image in degree
787 :param dec0: declination at the center of the image in degree
789 :return psf_r50: PSF r50 (flux radius 50%) evaluated at every object's position
790 :return psf_beta: PSF beta parameter evaluated at every object's position
791 :return psf_e1: PSF ellipticity 1-component evaluated at every object's position
792 :return psf_e2: PSF ellipticity 2-component evaluated at every object's position
793 """
795 def __str__(self):
796 return NAME
798 def __call__(self):
799 generator = PSF_GENERATOR[self.ctx.parameters.psf_type]
800 LOGGER.info(f"Generating PSF with {self.ctx.parameters.psf_type}")
802 psf_fwhm = np.array([])
804 def add_psf_to_objects(obj, ctx, name_str):
805 par = ctx.parameters
807 generator(name_str, ctx)
809 for p in ctx.psf_column_names:
810 # set catalog precision
811 arr = getattr(obj, p).astype(par.catalog_precision)
812 setattr(obj, p, arr)
814 if hasattr(self.ctx, "galaxies"):
815 add_psf_to_objects(self.ctx.galaxies, self.ctx, "galaxies")
816 update_psf_for_current_filter(self.ctx.galaxies, self.ctx)
818 self.ctx.galaxies.psf_fwhm = np.random.normal(
819 self.ctx.galaxies.psf_fwhm, self.ctx.parameters.psf_fwhm_variation_sigma
820 )
821 select = self.ctx.galaxies.psf_fwhm < 0
822 if np.any(select):
823 LOGGER.warning("Negative PSF FWHM detected. Setting to 0.")
824 self.ctx.galaxies.psf_fwhm[select] = 0.0
826 psf_fwhm = np.append(psf_fwhm, self.ctx.galaxies.psf_fwhm)
828 if hasattr(self.ctx, "stars"):
829 add_psf_to_objects(self.ctx.stars, self.ctx, "stars")
830 update_psf_for_current_filter(self.ctx.stars, self.ctx)
832 self.ctx.stars.psf_fwhm = np.random.normal(
833 self.ctx.stars.psf_fwhm, self.ctx.parameters.psf_fwhm_variation_sigma
834 )
835 select = self.ctx.stars.psf_fwhm < 0
836 if np.any(select):
837 LOGGER.warning("Negative PSF FWHM detected. Setting to 0.")
838 self.ctx.stars.psf_fwhm[select] = 0.0
839 psf_fwhm = np.append(psf_fwhm, self.ctx.stars.psf_fwhm)
841 if psf_fwhm.size > 0:
842 self.ctx.average_seeing = np.mean(psf_fwhm)
843 if not np.isfinite(self.ctx.average_seeing): # pragma: no cover
844 raise ValueError("Average seeing is not finite")
846 try:
847 del self.ctx.psf_r50_map
848 del self.ctx.psf_e1_map
849 del self.ctx.psf_e2_map
850 except AttributeError:
851 pass