From 9806d8ec1c44b52ccbc152599c2bc0de1f8c8b80 Mon Sep 17 00:00:00 2001 From: = Date: Fri, 14 Dec 2018 13:46:33 -0300 Subject: [PATCH] Change from __init__.py to own module. Minor optimization in logo reading in workers. Fix logo position. Adding logging. --- pcigale_plots/parsers/sed.py | 277 +++++++++++++++++++++++++++++++++++ 1 file changed, 277 insertions(+) create mode 100644 pcigale_plots/parsers/sed.py diff --git a/pcigale_plots/parsers/sed.py b/pcigale_plots/parsers/sed.py new file mode 100644 index 00000000..030f039b --- /dev/null +++ b/pcigale_plots/parsers/sed.py @@ -0,0 +1,277 @@ + +import logging +from itertools import repeat +from collections import OrderedDict + +from astropy.table import Table +import matplotlib + +matplotlib.use('Agg') +import matplotlib.pyplot as plt +import multiprocessing as mp +import numpy as np +import os +import pkg_resources +from scipy.constants import c +from pcigale.data import Database +from pcigale.utils import read_table +import matplotlib.gridspec as gridspec + +logger = logging.getLogger(__name__) + +# Wavelength limits (restframe) when plotting the best SED. +PLOT_L_MIN = 0.1 +PLOT_L_MAX = 5e5 + + +def sed(config, sed_type, best_results_file, nologo): + """Plot the best SED with associated observed and modelled fluxes. + """ + obs = read_table(config.configuration['data_file']) + mod = Table.read(best_results_file) + + with Database() as base: + filters = OrderedDict([(name, base.get_filter(name)) + for name in config.configuration['bands'] + if not (name.endswith('_err') or name.startswith('line')) ]) + + logo = False if nologo else plt.imread(pkg_resources.resource_filename(__name__, + "../resources/CIGALE.png")) + + with mp.Pool(processes=config.configuration['cores']) as pool: + pool.starmap(_sed_worker, zip(obs, mod, repeat(filters), + repeat(sed_type), repeat(logo))) + pool.close() + pool.join() + + +def _sed_worker(obs, mod, filters, sed_type, logo): + """Plot the best SED with the associated fluxes in bands + + Parameters + ---------- + obs: Table row + Data from the input file regarding one object. + mod: Table row + Data from the best model of one object. + filters: ordered dictionary of Filter objects + The observed fluxes in each filter. + sed_type: string + Type of SED to plot. It can either be "mJy" (flux in mJy and observed + frame) or "lum" (luminosity in W and rest frame) + logo: numpy.array | boolean + The readed logo image data. Has shape + (M, N) for grayscale images. + (M, N, 3) for RGB images. + (M, N, 4) for RGBA images. + Do not add the logo when set to False. + + """ + logger.debug("Starting worker") + + if os.path.isfile("out/{}_best_model.fits".format(obs['id'])): + + sed = Table.read("out/{}_best_model.fits".format(obs['id'])) + + filters_wl = np.array([filt.pivot_wavelength + for filt in filters.values()]) + wavelength_spec = sed['wavelength'] + obs_fluxes = np.array([obs[filt] for filt in filters.keys()]) + obs_fluxes_err = np.array([obs[filt+'_err'] + for filt in filters.keys()]) + mod_fluxes = np.array([mod["best."+filt] for filt in filters.keys()]) + if obs['redshift'] >= 0: + z = obs['redshift'] + else: # Redshift mode + z = mod['best.universe.redshift'] + DL = mod['best.universe.luminosity_distance'] + + if sed_type == 'lum': + xmin = PLOT_L_MIN + xmax = PLOT_L_MAX + + k_corr_SED = 1e-29 * (4.*np.pi*DL*DL) * c / (filters_wl*1e-9) + obs_fluxes *= k_corr_SED + obs_fluxes_err *= k_corr_SED + mod_fluxes *= k_corr_SED + + for cname in sed.colnames[1:]: + sed[cname] *= wavelength_spec + + filters_wl /= 1. + z + wavelength_spec /= 1. + z + elif sed_type == 'mJy': + xmin = PLOT_L_MIN * (1. + z) + xmax = PLOT_L_MAX * (1. + z) + k_corr_SED = 1. + + for cname in sed.colnames[1:]: + sed[cname] *= (wavelength_spec * 1e29 / + (c / (wavelength_spec * 1e-9)) / + (4. * np.pi * DL * DL)) + else: + logger.error("Unknown plot type") + + filters_wl /= 1000. + wavelength_spec /= 1000. + + wsed = np.where((wavelength_spec > xmin) & (wavelength_spec < xmax)) + figure = plt.figure() + gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1]) + if (sed.columns[1][wsed] > 0.).any(): + ax1 = plt.subplot(gs[0]) + ax2 = plt.subplot(gs[1]) + + # Stellar emission + if 'nebular.absorption_young' in sed.columns: + ax1.loglog(wavelength_spec[wsed], + (sed['stellar.young'][wsed] + + sed['attenuation.stellar.young'][wsed] + + sed['nebular.absorption_young'][wsed] + + sed['stellar.old'][wsed] + + sed['attenuation.stellar.old'][wsed] + + sed['nebular.absorption_old'][wsed]), + label="Stellar attenuated ", color='orange', + marker=None, nonposy='clip', linestyle='-', + linewidth=0.5) + else: + ax1.loglog(wavelength_spec[wsed], + (sed['stellar.young'][wsed] + + sed['attenuation.stellar.young'][wsed] + + sed['stellar.old'][wsed] + + sed['attenuation.stellar.old'][wsed]), + label="Stellar attenuated ", color='orange', + marker=None, nonposy='clip', linestyle='-', + linewidth=0.5) + ax1.loglog(wavelength_spec[wsed], + (sed['stellar.old'][wsed] + + sed['stellar.young'][wsed]), + label="Stellar unattenuated", color='b', marker=None, + nonposy='clip', linestyle='--', linewidth=0.5) + # Nebular emission + if 'nebular.lines_young' in sed.columns: + ax1.loglog(wavelength_spec[wsed], + (sed['nebular.lines_young'][wsed] + + sed['nebular.lines_old'][wsed] + + sed['nebular.continuum_young'][wsed] + + sed['nebular.continuum_old'][wsed] + + sed['attenuation.nebular.lines_young'][wsed] + + sed['attenuation.nebular.lines_old'][wsed] + + sed['attenuation.nebular.continuum_young'][wsed] + + sed['attenuation.nebular.continuum_old'][wsed]), + label="Nebular emission", color='y', marker=None, + nonposy='clip', linewidth=.5) + # Dust emission Draine & Li + if 'dust.Umin_Umin' in sed.columns: + ax1.loglog(wavelength_spec[wsed], + (sed['dust.Umin_Umin'][wsed] + + sed['dust.Umin_Umax'][wsed]), + label="Dust emission", color='r', marker=None, + nonposy='clip', linestyle='-', linewidth=0.5) + # Dust emission Dale + if 'dust' in sed.columns: + ax1.loglog(wavelength_spec[wsed], sed['dust'][wsed], + label="Dust emission", color='r', marker=None, + nonposy='clip', linestyle='-', linewidth=0.5) + # AGN emission Fritz + if 'agn.fritz2006_therm' in sed.columns: + ax1.loglog(wavelength_spec[wsed], + (sed['agn.fritz2006_therm'][wsed] + + sed['agn.fritz2006_scatt'][wsed] + + sed['agn.fritz2006_agn'][wsed]), + label="AGN emission", color='g', marker=None, + nonposy='clip', linestyle='-', linewidth=0.5) + # Radio emission + if 'radio_nonthermal' in sed.columns: + ax1.loglog(wavelength_spec[wsed], + sed['radio_nonthermal'][wsed], + label="Radio nonthermal", color='brown', + marker=None, nonposy='clip', linestyle='-', + linewidth=0.5) + + ax1.loglog(wavelength_spec[wsed], sed['L_lambda_total'][wsed], + label="Model spectrum", color='k', nonposy='clip', + linestyle='-', linewidth=1.5) + + ax1.set_autoscale_on(False) + s = np.argsort(filters_wl) + filters_wl = filters_wl[s] + mod_fluxes = mod_fluxes[s] + obs_fluxes = obs_fluxes[s] + obs_fluxes_err = obs_fluxes_err[s] + ax1.scatter(filters_wl, mod_fluxes, marker='o', color='r', s=8, + zorder=3, label="Model fluxes") + mask_ok = np.logical_and(obs_fluxes > 0., obs_fluxes_err > 0.) + ax1.errorbar(filters_wl[mask_ok], obs_fluxes[mask_ok], + yerr=obs_fluxes_err[mask_ok]*3, ls='', marker='s', + label='Observed fluxes', markerfacecolor='None', + markersize=6, markeredgecolor='b', capsize=0.) + mask_uplim = np.logical_and(np.logical_and(obs_fluxes > 0., + obs_fluxes_err < 0.), + obs_fluxes_err > -9990. * k_corr_SED) + if not mask_uplim.any() == False: + ax1.errorbar(filters_wl[mask_uplim], obs_fluxes[mask_uplim], + yerr=obs_fluxes_err[mask_uplim]*3, ls='', + marker='v', label='Observed upper limits', + markerfacecolor='None', markersize=6, + markeredgecolor='g', capsize=0.) + mask_noerr = np.logical_and(obs_fluxes > 0., + obs_fluxes_err < -9990. * k_corr_SED) + if not mask_noerr.any() == False: + ax1.errorbar(filters_wl[mask_noerr], obs_fluxes[mask_noerr], + ls='', marker='s', markerfacecolor='None', + markersize=6, markeredgecolor='r', + label='Observed fluxes, no errors', capsize=0.) + mask = np.where(obs_fluxes > 0.) + ax2.errorbar(filters_wl[mask], + (obs_fluxes[mask]-mod_fluxes[mask])/obs_fluxes[mask], + yerr=obs_fluxes_err[mask]/obs_fluxes[mask]*3, + marker='_', label="(Obs-Mod)/Obs", color='k', + capsize=0.) + ax2.plot([xmin, xmax], [0., 0.], ls='--', color='k') + ax2.set_xscale('log') + ax2.minorticks_on() + + figure.subplots_adjust(hspace=0., wspace=0.) + + ax1.set_xlim(xmin, xmax) + ymin = min(np.min(obs_fluxes[mask_ok]), + np.min(mod_fluxes[mask_ok])) + if not mask_uplim.any() == False: + ymax = max(max(np.max(obs_fluxes[mask_ok]), + np.max(obs_fluxes[mask_uplim])), + max(np.max(mod_fluxes[mask_ok]), + np.max(mod_fluxes[mask_uplim]))) + else: + ymax = max(np.max(obs_fluxes[mask_ok]), + np.max(mod_fluxes[mask_ok])) + ax1.set_ylim(1e-1*ymin, 1e1*ymax) + ax2.set_xlim(xmin, xmax) + ax2.set_ylim(-1.0, 1.0) + if sed_type == 'lum': + ax2.set_xlabel("Rest-frame wavelength [$\mu$m]") + ax1.set_ylabel("Luminosity [W]") + ax2.set_ylabel("Relative residual luminosity") + else: + ax2.set_xlabel("Observed wavelength [$\mu$m]") + ax1.set_ylabel("Flux [mJy]") + ax2.set_ylabel("Relative residual flux") + ax1.legend(fontsize=6, loc='best', fancybox=True, framealpha=0.5) + ax2.legend(fontsize=6, loc='best', fancybox=True, framealpha=0.5) + plt.setp(ax1.get_xticklabels(), visible=False) + plt.setp(ax1.get_yticklabels()[1], visible=False) + figure.suptitle("Best model for {} at z = {}. Reduced $\chi^2$={}". + format(obs['id'], np.round(z, decimals=3), + np.round(mod['best.reduced_chi_square'], + decimals=2))) + if logo is not False: + figure_height = figure.get_figheight() * figure.dpi + figure.figimage(logo, 12, figure_height - 67, origin='upper', zorder=0, + alpha=1) + + figure.savefig("out/{}_best_model.pdf".format(obs['id'])) + plt.close(figure) + else: + logger.error("No valid best SED found for {}. No plot created.".format(obs['id'])) + else: + logger.error("No SED found for {}. No plot created.".format(obs['id'])) -- GitLab