From ac8ad5516ab8fa642080370d0c6feb8de511d3ac Mon Sep 17 00:00:00 2001 From: Robel Geda Date: Wed, 11 Jan 2023 15:00:38 -0500 Subject: [PATCH] plot update --- petrofit/petrosian/core.py | 261 +++++++++++++++++++++++-------- petrofit/petrosian/correction.py | 16 +- 2 files changed, 206 insertions(+), 71 deletions(-) diff --git a/petrofit/petrosian/core.py b/petrofit/petrosian/core.py index 5c5e67f..c6502df 100644 --- a/petrofit/petrosian/core.py +++ b/petrofit/petrosian/core.py @@ -344,28 +344,53 @@ def fraction_to_r(fraction, r_list, flux_list, r_petrosian, return r_frac, r_frac_err -def calculate_concentration_index(r_list, flux_list, r_total_flux, fraction_1=0.2, fraction_2=0.8): +def calculate_concentration_index(fraction_1, fraction_2, + r_list, flux_list, r_petrosian, + flux_err=None, r_petrosian_err=None, + epsilon=2., epsilon_fraction=0.99, + interp_kind='cubic', interp_num=5000): """ Calculates Petrosian concentration index. - ``concentration_index = 5 * np.log10( r(fraction_2) / r(fraction_1) )`` + ``concentration_index = C_1_2 = 5 * np.log10( r(fraction_2) / r(fraction_1) )`` Parameters ---------- + fraction_1 : float + Fraction of total light enclosed by the radius in the denominator. + + fraction_2 : float + Fraction of total light enclosed by the radius in the numerator. + r_list : numpy.array Array of radii in pixels. flux_list : numpy.array Array of photometric flux values. - r_total_flux : float - Total flux radius (hint: `petrofit.petrosian.calculate_r_total_flux` can be used to measure this). + r_petrosian : float + Petrosian radius - fraction_1 : float - Fraction of total light enclosed by the radius in the denominator. + flux_err : int or numpy.array, optional + Array of errors in the flux values. If None, erros are not computed even if + `r_petro_err` is provided. - fraction_2 : float - Fraction of total light enclosed by the radius in the numerator. + r_petrosian_err : float, optional + 1-sigma error in r_petro. If set to `None` and `flux_err` is provided, + `r_petro_err` is assumed to be 0 when computing errors. + + epsilon : float, default=2. + Epsilon value (used to determine `r_epsilon`). + N.B: `r_epsilon = r_petrosian` * epsilon` + + epsilon_fraction: float, default=0.99 + Fraction of total flux that is recovered by `r_petrosian * epsilon`. + + interp_kind : str or int, optional + Specifies the kind of interpolation used on the curve of growth (i.e `r_list` vs `flux_list`) + + interp_num : int + Number of interpolation function sampling radii. Returns ------- @@ -381,14 +406,21 @@ def calculate_concentration_index(r_list, flux_list, r_total_flux, fraction_1=0. Concentration index """ - if r_total_flux > max(r_list): - return [np.nan, np.nan, np.nan] + r1 = fraction_to_r(fraction_1, r_list, flux_list, r_petrosian, + flux_err=flux_err, r_petrosian_err=r_petrosian_err, + epsilon=epsilon, epsilon_fraction=epsilon_fraction, + interp_kind=interp_kind, interp_num=interp_num) + + r2 = fraction_to_r(fraction_2, r_list, flux_list, r_petrosian, + flux_err=flux_err, r_petrosian_err=r_petrosian_err, + epsilon=epsilon, epsilon_fraction=epsilon_fraction, + interp_kind=interp_kind, interp_num=interp_num) - r1 = fraction_flux_to_r(r_list, flux_list, r_total_flux, fraction=fraction_1) - r2 = fraction_flux_to_r(r_list, flux_list, r_total_flux, fraction=fraction_2) + r1 = r1[0] + r2 = r2[0] - # if np.any(np.isnan(np.array([r1, r2]))): - # return [np.nan, np.nan, np.nan] + if np.any(np.isnan(np.array([r1, r2]))): + return [np.nan, np.nan, np.nan] return r1, r2, 5 * np.log10(r2 / r1) @@ -492,6 +524,52 @@ def _update_petrosian(self): flux_err=self.flux_err) self.has_petrosian_err = self.petrosian_err is not None + def _calculate_petrosian_r(self): + return calculate_petrosian_r(self.r_list, self.petrosian_list, + petrosian_err=self.petrosian_err, + eta=self.eta, + interp_kind=self.interp_kind, + interp_num=self.interp_num) + + def _calculate_fraction_to_r(self, fraction): + return fraction_to_r(fraction, self.r_list, self.flux_list, self.r_petrosian, + flux_err=self.flux_err, r_petrosian_err=self.r_petrosian_err, + epsilon=self.epsilon, epsilon_fraction=self.epsilon_fraction, + interp_kind=self.interp_kind, interp_num=self.interp_num) + + def concentration_index(self, fraction_1=0.2, fraction_2=0.8): + """ + Calculates Petrosian concentration index. + + ``concentration_index = C_1_2 = 5 * np.log10( r(fraction_2) / r(fraction_1) )`` + + Parameters + ---------- + fraction_1 : float + Fraction of total light enclosed by the radius in the denominator. + + fraction_2 : float + Fraction of total light enclosed by the radius in the numerator. + + Returns + ------- + r_fraction_1, r_fraction_2, concentration_index + + * r_fraction_1 : float or np.nan + Radius containing `fraction_1` of the total flux. + + * r_fraction_2: float or np.nan + Radius containing `fraction_2` of the total flux. + + * concentration_index : float or np.nan + Concentration index + """ + return calculate_concentration_index(fraction_1, fraction_2, + self.r_list, self.flux_list, self.r_petrosian, + flux_err=self.flux_err, r_petrosian_err=self.r_petrosian_err, + epsilon=self.epsilon, epsilon_fraction=self.epsilon_fraction, + interp_kind=self.interp_kind, interp_num=self.interp_num) + @property def r_list(self): return self._r_list @@ -548,19 +626,6 @@ def total_flux_fraction(self): def total_flux_fraction(self, value): self._total_flux_fraction = value - def _calculate_petrosian_r(self): - return calculate_petrosian_r(self.r_list, self.petrosian_list, - petrosian_err=self.petrosian_err, - eta=self.eta, - interp_kind=self.interp_kind, - interp_num=self.interp_num) - - def _calculate_fraction_to_r(self, fraction): - return fraction_to_r(fraction, self.r_list, self.flux_list, self.r_petrosian, - flux_err=self.flux_err, r_petrosian_err=self.r_petrosian_err, - epsilon=self.epsilon, epsilon_fraction=self.epsilon_fraction, - interp_kind=self.interp_kind) - @property def r_petrosian(self): """ @@ -655,36 +720,6 @@ def r_total_flux_arcsec(self, wcs): return pixel_to_angular(r_total_flux, wcs).value return np.nan - def concentration_index(self, fraction_1=0.2, fraction_2=0.8): - """ - Calculates Petrosian concentration index. - - ``concentration_index = 5 * np.log10( r(fraction_2) / r(fraction_1) )`` - - Parameters - ---------- - fraction_1 : float - Fraction of total light enclosed by the radius in the denominator. - - fraction_2 : float - Fraction of total light enclosed by the radius in the numerator. - - Returns - ------- - r_fraction_1, r_fraction_2, concentration_index - - * r_fraction_1 : float or np.nan - Radius containing `fraction_1` of the total flux. - - * r_fraction_2: float or np.nan - Radius containing `fraction_2` of the total flux. - - * concentration_index : float or np.nan - Concentration index - """ - return calculate_concentration_index(self.r_list, self.flux_list, self.r_total_flux, - fraction_1=fraction_1, fraction_2=fraction_2) - @property def c2080(self): """``c2080 = 5 * np.log10(r_80 / r_20)``""" @@ -695,7 +730,34 @@ def c5090(self): """``c5090 = 5 * np.log10(r_90 / r_50)``""" return self.concentration_index(fraction_1=0.5, fraction_2=0.9)[-1] - def plot(self, plot_r=False, plot_normalized_flux=False): + def _plot_radii(self, ax, radius_unit='pix', err_capsize=3, ): + radius_unit = '' if radius_unit is None else str(radius_unit) + + r_petrosian = self.r_petrosian + if not np.isnan(r_petrosian): + ax.axvline(r_petrosian, linestyle='--', + label="$r(\eta_{{{}}})={:0.4f}$ {}".format(self.eta, r_petrosian, radius_unit)) + + r_epsilon = self.r_petrosian * self.epsilon + if not np.isnan(r_epsilon): + epsilon_fraction = int(self.epsilon_fraction * 100) + ax.axvline(r_epsilon, linestyle='--', c='green', + label="$r_{{\epsilon}}(L_{{{}}}, \epsilon = {}) = {:0.4f}$ {}".format(epsilon_fraction, + self.epsilon, r_epsilon, + radius_unit)) + + r_total_flux = self.r_total_flux + if not np.isnan(r_total_flux): + total_flux_fraction = int(self.total_flux_fraction * 100) + ax.axvline(r_total_flux, linestyle='--', c='black', + label="$r_{{total}}(L_{{{}}}) = {:0.4f}$ {}".format(total_flux_fraction, r_total_flux, + radius_unit)) + + def plot(self, plot_r=False, show_normalized_flux=False, + title='Petrosian Profile', radius_unit='pix', + ax=None, err_alpha=0.2, err_capsize=3, + show_legend=True, legend_fontsize=None, + ax_fontsize=None, tick_fontsize=None): """ Plots Petrosian profile. @@ -707,10 +769,83 @@ def plot(self, plot_r=False, plot_normalized_flux=False): plot_normalized_flux: Over-plot the flux curve of growth by normalizing the flux axis (max_flux=1). """ - plot_petrosian(self.r_list, self.area_list, self.flux_list, epsilon=self.epsilon, eta=self.eta, plot_r=plot_r) + radius_unit = '' if radius_unit is None else str(radius_unit) + if ax is None: + ax = plt.gca() + + ax.errorbar(self.r_list, self.petrosian_list, yerr=self.petrosian_err, + marker='o', capsize=err_capsize, + label="Data") + + if err_alpha is not None and self.has_petrosian_err and err_alpha > 0: + ax.fill_between(self.r_list, + self.petrosian_list - self.petrosian_err, + self.petrosian_list + self.petrosian_err, + alpha=err_alpha) + + r_petrosian = self.r_petrosian + r_petrosian_err = self.r_petrosian_err + if not np.isnan(r_petrosian) and not np.isnan(r_petrosian_err): + ax.errorbar(r_petrosian, self.eta, xerr=r_petrosian_err, + marker='o', capsize=err_capsize) + ax.axhline(self.eta, linestyle='--') + + ax.axhline(0, c='black') + + self._plot_radii(ax, radius_unit=radius_unit, err_capsize=err_capsize) + + ax.set_title(title, fontsize=ax_fontsize) + ax.set_xlabel("Aperture Radius" + " [{}]".format(radius_unit) if radius_unit else "", fontsize=ax_fontsize) + ax.set_ylabel("Petrosian Value $\eta(r)$", fontsize=ax_fontsize) + + ax.tick_params(axis='both', which='major', labelsize=tick_fontsize) + ax.tick_params(axis='both', which='minor', labelsize=tick_fontsize) + + if show_legend: + ax.legend(fontsize=legend_fontsize) + + return ax + + def plot_cog(self, plot_r=False, show_normalized_flux=False, + title='Petrosian Profile', radius_unit='pix', + ax=None, err_alpha=0.2, err_capsize=3, + show_legend=True, legend_fontsize=None, + ax_fontsize=None, tick_fontsize=None): + """ + """ + radius_unit = '' if radius_unit is None else str(radius_unit) + if ax is None: + ax = plt.gca() + + ax.errorbar(self.r_list, self.flux_list, yerr=self.flux_err, + marker='o', capsize=err_capsize, + label="Data") - if plot_normalized_flux: - plt.plot(self.r_list, self.flux_list / self.flux_list.max(), label='Normalized Flux', linestyle='--') + if err_alpha is not None and self.has_petrosian_err and err_alpha > 0: + ax.fill_between(self.r_list, + self.flux_list - self.flux_err, + self.flux_list + self.flux_err, + alpha=err_alpha) + + self._plot_radii(ax, radius_unit=radius_unit, err_capsize=err_capsize) + + total_flux = self.total_flux + if not np.isnan(total_flux): + ax.axhline(total_flux, linestyle='--', c='black', ) + + ax.axhline(0, c='black') + + ax.set_title(title, fontsize=ax_fontsize) + ax.set_xlabel("Aperture Radius" + " [{}]".format(radius_unit) if radius_unit else "", fontsize=ax_fontsize) + ax.set_ylabel("Petrosian Value $\eta(r)$", fontsize=ax_fontsize) + + ax.tick_params(axis='both', which='major', labelsize=tick_fontsize) + ax.tick_params(axis='both', which='minor', labelsize=tick_fontsize) + + if show_legend: + ax.legend(fontsize=legend_fontsize) + + return ax def imshow(self, position=(0, 0), elong=1., theta=0., color=None, lw=None): """ @@ -735,7 +870,8 @@ def imshow(self, position=(0, 0), elong=1., theta=0., color=None, lw=None): """ labels = ['r_half_light', 'r_total_flux', 'r_20', 'r_80'] - radii = [self.r_half_light, self.r_total_flux, self.fraction_flux_to_r(.2), self.fraction_flux_to_r(.8)] + radii = [self.r_half_light, self.r_total_flux, self._calculate_fraction_to_r(.2)[0], + self._calculate_fraction_to_r(.8)[0]] colors = ['r', 'r', 'b', 'b'] linestyles = ['dashed', 'solid', 'dotted', 'dashdot'] @@ -745,4 +881,3 @@ def imshow(self, position=(0, 0), elong=1., theta=0., color=None, lw=None): label=label, linestyle=ls, color=color if color else default_color, lw=lw) plt.scatter(*position, marker='+', color=color if color else 'red') - diff --git a/petrofit/petrosian/correction.py b/petrofit/petrosian/correction.py index 3e11a61..41f2158 100644 --- a/petrofit/petrosian/correction.py +++ b/petrofit/petrosian/correction.py @@ -189,12 +189,12 @@ def _generate_petrosian_correction(args): amplitude = 100 / np.exp(gammaincinv(2. * n, 0.5)) # Total flux - r_100 = sersic_enclosed( + L_total = sersic_enclosed( np.inf, amplitude=amplitude, r_eff=r_eff, n=n) - total_flux = r_100 * 0.99 + total_flux = L_total * 0.99 # Calculate radii r_20, r_80, r_total_flux = [sersic_enclosed_inv( @@ -259,7 +259,7 @@ def _generate_petrosian_correction(args): # Compute corrections corrected_epsilon = model_r_total_flux / p.r_petrosian - corrected_epsilon80 = model_r_80 / p.r_petrosian + corrected_epsilon_80 = model_r_80 / p.r_petrosian corrected_p = copy(p) corrected_p.epsilon = corrected_epsilon @@ -267,8 +267,8 @@ def _generate_petrosian_correction(args): # Make output list # ---------------- # Petrosian indices - petrosian_list = calculate_petrosian(p.area_list, p.flux_list) - p02, p03, p04, p05 = [calculate_petrosian_r(p.r_list, petrosian_list, eta=i) for i in (0.2, 0.3, 0.4, 0.5)] + petrosian_list = calculate_petrosian(p.area_list, p.flux_list)[0] + p02, p03, p04, p05 = [calculate_petrosian_r(p.r_list, petrosian_list, eta=i)[0] for i in (0.2, 0.3, 0.4, 0.5)] assert np.round(p.r_petrosian, 6) == np.round(p02, 6) u_r_eff = p.fraction_flux_to_r(0.5) @@ -279,10 +279,10 @@ def _generate_petrosian_correction(args): c_r_20 = corrected_p.fraction_flux_to_r(0.2) c_r_80 = corrected_p.fraction_flux_to_r(0.8) - row = [n, r_eff, r_20, r_80, r_total_flux, r_100, + row = [n, r_eff, r_20, r_80, r_total_flux, L_total, p02, p03, p04, p05, 5 * np.log(p02 / p05), 5 * np.log(p02 / p03), p.epsilon, u_r_80 / p.r_petrosian, u_r_eff, p.r_total_flux, u_r_20, u_r_80, p.c2080, - corrected_epsilon, corrected_epsilon80, c_r_eff, corrected_p.r_total_flux, c_r_20, c_r_80, + corrected_epsilon, corrected_epsilon_80, c_r_eff, corrected_p.r_total_flux, c_r_20, c_r_80, corrected_p.c2080, ] if plot: corrected_p.plot(True, True) @@ -382,7 +382,7 @@ def generate_petrosian_sersic_correction(output_yaml_name, psf=None, r_eff_list= rows = ProgressBar.map(_generate_petrosian_correction, args, multiprocess=n_cpu, ipython_widget=ipython_widget, step=step) - names = ['n', 'r_eff', 'sersic_r_20', 'sersic_r_80', 'sersic_r_99', 'sersic_r_100', + names = ['n', 'r_eff', 'sersic_r_20', 'sersic_r_80', 'sersic_r_99', 'sersic_L_inf', 'p02', 'p03', 'p04', 'p05', 'p0502', 'p0302', 'u_epsilon', 'u_epsilon_80', 'u_r_eff', 'u_r_99', 'u_r_20', 'u_r_80', 'u_c2080', 'c_epsilon', 'c_epsilon_80', 'c_r_eff', 'c_r_99', 'c_r_20', 'c_r_80', 'c_c2080', ]