From dd7566d9b6b9a3697efb34cf362cd01634987cd7 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Thu, 30 Apr 2020 10:17:54 +0200 Subject: [PATCH] Object support for duplicate filter names --- species/analysis/photometry.py | 47 +++------ species/core/box.py | 21 ++-- species/data/companions.py | 6 +- species/data/database.py | 151 +++++++++++++++++++++-------- species/plot/plot_spectrum.py | 136 +++++++++++++++++--------- species/util/phot_util.py | 28 ++++-- test/test_read/test_calibration.py | 4 +- test/test_read/test_isochrone.py | 8 +- test/test_read/test_model.py | 20 ++-- test/test_read/test_planck.py | 8 +- 10 files changed, 270 insertions(+), 159 deletions(-) diff --git a/species/analysis/photometry.py b/species/analysis/photometry.py index 8373e124..2a6ef319 100644 --- a/species/analysis/photometry.py +++ b/species/analysis/photometry.py @@ -131,19 +131,6 @@ def spectrum_to_flux(self, raise ValueError('Calculation of the mean flux is not possible because the ' 'wavelength array is empty.') - indices = np.where((wavelength > self.wavel_range[0]) & - (wavelength < self.wavel_range[1]))[0] - - if indices.size == 1: - raise ValueError('Calculating synthetic photometry requires more than one ' - 'wavelength point.') - - wavelength = wavelength[indices] - flux = flux[indices] - - if error is not None: - error = error[indices] - indices = np.where((self.wavel_range[0] <= wavelength) & (wavelength <= self.wavel_range[1]))[0] @@ -160,7 +147,7 @@ def spectrum_to_flux(self, warnings.warn(f'The filter profile of {self.filter_name} ' f'({self.wavel_range[0]:.4f}-{self.wavel_range[1]:.4f}) extends ' - f' beyond the wavelength range of the spectrum ({wavelength[0]:.4f} ' + f'beyond the wavelength range of the spectrum ({wavelength[0]:.4f} ' f'-{wavelength[-1]:.4f}). The flux is set to NaN. Setting the ' f'\'threshold\' parameter will loosen the wavelength constraints.') @@ -202,21 +189,19 @@ def spectrum_to_flux(self, syn_flux = integral1/integral2 if error is not None and not np.any(np.isnan(error)): - error_flux = np.zeros(200) + phot_random = np.zeros(200) for i in range(200): - spec_random = flux+np.random.normal(loc=0., - scale=1., - size=wavelength.shape[0])*error - - spec_tmp = self.spectrum_to_flux(wavelength, - spec_random, - error=None, - threshold=threshold)[0] + spec_random = flux + np.random.normal(loc=0., + scale=1., + size=wavelength.shape[0])*error - error_flux[i] = spec_tmp + phot_random[i] = self.spectrum_to_flux(wavelength, + spec_random, + error=None, + threshold=threshold)[0] - error_flux = np.std(error_flux) + error_flux = np.std(phot_random) else: error_flux = None @@ -269,21 +254,21 @@ def spectrum_to_magnitude(self, app_mag = self.vega_mag - 2.5*math.log10(syn_flux[0]/zp_flux) if error is not None and not np.any(np.isnan(error)): - error_app_mag = np.zeros(200) + mag_random = np.zeros(200) for i in range(200): - spec_random = flux+np.random.normal(loc=0., - scale=1., - size=wavelength.shape[0])*error + spec_random = flux + np.random.normal(loc=0., + scale=1., + size=wavelength.shape[0])*error flux_random = self.spectrum_to_flux(wavelength, spec_random, error=None, threshold=threshold) - error_app_mag[i] = self.vega_mag - 2.5*np.log10(flux_random[0]/zp_flux) + mag_random[i] = self.vega_mag - 2.5*np.log10(flux_random[0]/zp_flux) - error_app_mag = np.std(error_app_mag) + error_app_mag = np.std(mag_random) else: error_app_mag = None diff --git a/species/core/box.py b/species/core/box.py index 69debbd6..b27a7765 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -60,13 +60,20 @@ def create_box(boxtype, elif boxtype == 'photometry': box = PhotometryBox() - box.name = kwargs['name'] - box.sptype = kwargs['sptype'] - box.wavelength = kwargs['wavelength'] - box.flux = kwargs['flux'] - box.app_mag = kwargs['app_mag'] - box.abs_mag = kwargs['abs_mag'] - box.filter_name = kwargs['filter_name'] + if 'name' in kwargs: + box.name = kwargs['name'] + if 'sptype' in kwargs: + box.sptype = kwargs['sptype'] + if 'wavelength' in kwargs: + box.wavelength = kwargs['wavelength'] + if 'flux' in kwargs: + box.flux = kwargs['flux'] + if 'app_mag' in kwargs: + box.app_mag = kwargs['app_mag'] + if 'abs_mag' in kwargs: + box.abs_mag = kwargs['abs_mag'] + if 'filter_name' in kwargs: + box.filter_name = kwargs['filter_name'] elif boxtype == 'residuals': box = ResidualsBox() diff --git a/species/data/companions.py b/species/data/companions.py index 036bdc5e..ad3bd3af 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -102,10 +102,8 @@ def get_data(): 'PDS 70 b': {'distance': (113.43, 0.52), 'app_mag': {'Paranal/SPHERE.IRDIS_D_H23_2': (17.94, 0.24), # Keppler et al. 2018 'Paranal/SPHERE.IRDIS_D_H23_3': (17.95, 0.17), # Keppler et al. 2018 - 'Paranal/SPHERE.IRDIS_D_K12_1': (16.78, 0.31), # Stolker et al. in prep. - 'Paranal/SPHERE.IRDIS_D_K12_2': (16.23, 0.32), # Stolker et al. in prep. - # 'Paranal/SPHERE.IRDIS_D_K12_1': (16.68, 0.04), # Stolker et al. in prep. - # 'Paranal/SPHERE.IRDIS_D_K12_2': (16.35, 0.07), # Stolker et al. in prep. + 'Paranal/SPHERE.IRDIS_D_K12_1': [(16.78, 0.31), (16.68, 0.04)], # Stolker et al. in prep. + 'Paranal/SPHERE.IRDIS_D_K12_2': [(16.23, 0.32), (16.35, 0.07)], # Stolker et al. in prep. 'Paranal/NACO.Lp': (14.08, 0.33), # Stolker et al. in prep. 'Paranal/NACO.NB405': (13.91, 0.34), # Stolker et al. in prep. 'Paranal/NACO.Mp': (13.64, 0.22)}}, # Stolker et al. in prep. diff --git a/species/data/database.py b/species/data/database.py index 7b60069f..0473c792 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -7,11 +7,14 @@ import warnings import configparser +from typing import Optional, Union, List, Tuple, Dict + import h5py import tqdm import emcee import numpy as np +from typeguard import typechecked from astropy.io import fits # from petitRADTRANS import Radtrans @@ -128,17 +131,19 @@ def delete_data(self, else: warnings.warn(f'The dataset {dataset} is not found in {self.database}.') + @typechecked def add_companion(self, - name=None): + name: Union[Optional[str], Optional[List[str]]]) -> None: """ Function for adding the magnitudes of directly imaged planets and brown dwarfs from :class:`~species.data.companions.get_data` to the database. Parameters ---------- - name : list(str, ), None - List with names of the directly imaged planets and brown dwarfs (e.g. ``['HR 8799 b', - '51 Eri b', 'PZ Tel B']``). All the available companion data are added if set to None. + name : str, list(str, ), None + Name or list with names of the directly imaged planets and brown dwarfs (e.g. + ``'HR 8799 b'`` or ``['HR 8799 b', '51 Eri b', 'PZ Tel B']``). All the available + companion data are added if set to ``None``. Returns ------- @@ -410,11 +415,17 @@ def add_model(self, h5_file.close() + @typechecked def add_object(self, - object_name, - distance=None, - app_mag=None, - spectrum=None): + object_name: str, + distance: Optional[Tuple[float, float]] = None, + app_mag: Optional[Dict[str, + Union[Tuple[float, float], + List[Tuple[float, float]]]]] = None, + spectrum: Optional[Dict[str, + Tuple[str, + Optional[str], + Optional[float]]]] = None) -> None: """ Function for adding the photometric and/or spectroscopic data of an object to the database. @@ -426,8 +437,10 @@ def add_object(self, Distance and uncertainty (pc). Not stored if set to None. app_mag : dict, None Dictionary with the filter names, apparent magnitudes, and uncertainties. For example, - ``{'Paranal/NACO.Lp': (15., 0.2), 'Paranal/NACO.Mp': (13., 0.3)}``. Not stored if set - to None. + ``{'Paranal/NACO.Lp': (15., 0.2), 'Paranal/NACO.Mp': (13., 0.3)}``. For the use of + duplicate filter names, the magnitudes have to be provided in a list, for example + ``{'Paranal/NACO.Lp': [(15., 0.2), (14.5, 0.5)], 'Paranal/NACO.Mp': (13., 0.3)}``. + No photometric data is stored if set to ``None``. spectrum : dict, None Dictionary with the spectrum, optional covariance matrix, and spectral resolution for each instrument. The input data can either have a FITS or ASCII format. The spectra @@ -475,37 +488,94 @@ def add_object(self, error = {} if app_mag is not None: - for item in app_mag: - try: - synphot = photometry.SyntheticPhotometry(item) - flux[item], error[item] = synphot.magnitude_to_flux(app_mag[item][0], - app_mag[item][1]) + for mag_item in app_mag: + if isinstance(app_mag[mag_item], tuple): - except KeyError: - warnings.warn(f'Filter \'{item}\' is not available on the SVO Filter Profile ' - f'Service so a flux calibration can not be done. Please add the ' - f'filter manually with the \'add_filter\' function. For now, ' - f'only the \'{item}\' magnitude of \'{object_name}\' is stored.') + try: + synphot = photometry.SyntheticPhotometry(mag_item) + flux[mag_item], error[mag_item] = synphot.magnitude_to_flux( + app_mag[mag_item][0], app_mag[mag_item][1]) - # Write NaNs if the filter is not available - flux[item], error[item] = np.nan, np.nan + except KeyError: + warnings.warn(f'Filter \'{mag_item}\' is not available on the SVO Filter ' + f'Profile Service so a flux calibration can not be done. Please ' + f'add the filter manually with the \'add_filter\' function. For ' + f'now, only the \'{mag_item}\' magnitude of \'{object_name}\' is ' + f'stored.') - for item in app_mag: - if f'objects/{object_name}/{item}' in h5_file: - del h5_file[f'objects/{object_name}/{item}'] + # Write NaNs if the filter is not available + flux[mag_item], error[mag_item] = np.nan, np.nan + + elif isinstance(app_mag[mag_item], list): + flux_list = [] + error_list = [] - print(f' - {item}:') - print(f' - Apparent magnitude = {app_mag[item][0]:.2f} +/- {app_mag[item][1]:.2f}') - print(f' - Flux (W m-2 um-1) = {flux[item]:.2e} +/- {error[item]:.2e}') + for i, dupl_item in enumerate(app_mag[mag_item]): - data = np.asarray([app_mag[item][0], - app_mag[item][1], - flux[item], - error[item]]) + try: + synphot = photometry.SyntheticPhotometry(mag_item) + flux_dupl, error_dupl = synphot.magnitude_to_flux( + dupl_item[0], dupl_item[1]) + + except KeyError: + warnings.warn(f'Filter \'{mag_item}\' is not available on the SVO Filter ' + f'Profile Service so a flux calibration can not be done. Please add the ' + f'filter manually with the \'add_filter\' function. For now, ' + f'only the \'{mag_item}\' magnitude of \'{object_name}\' is stored.') + + # Write NaNs if the filter is not available + flux_dupl, error_dupl = np.nan, np.nan + + flux_list.append(flux_dupl) + error_list.append(error_dupl) + + flux[mag_item] = flux_list + error[mag_item] = error_list + + else: + raise ValueError('The values in the dictionary with magnitudes should be ' + 'tuples or a list with tuples (in case duplicate filter ' + 'names are required).') + + for mag_item in app_mag: + if f'objects/{object_name}/{mag_item}' in h5_file: + del h5_file[f'objects/{object_name}/{mag_item}'] + + if isinstance(app_mag[mag_item], tuple): + n_phot = 1 + print(f' - {mag_item}:') + print(f' - Apparent magnitude = {app_mag[mag_item][0]:.2f} +/- {app_mag[mag_item][1]:.2f}') + print(f' - Flux (W m-2 um-1) = {flux[mag_item]:.2e} +/- {error[mag_item]:.2e}') + + data = np.asarray([app_mag[mag_item][0], + app_mag[mag_item][1], + flux[mag_item], + error[mag_item]]) + + elif isinstance(app_mag[mag_item], list): + n_phot = len(app_mag[mag_item]) + print(f' - {mag_item} ({n_phot} values):') + + mag_list = [] + mag_err_list = [] + + for i, dupl_item in enumerate(app_mag[mag_item]): + print(f' - Apparent magnitude = {app_mag[mag_item][i][0]:.2f} +/- {app_mag[mag_item][i][1]:.2f}') + print(f' - Flux (W m-2 um-1) = {flux[mag_item][i]:.2e} +/- {error[mag_item][i]:.2e}') + + mag_list.append(app_mag[mag_item][i][0]) + mag_err_list.append(app_mag[mag_item][i][1]) + + data = np.asarray([mag_list, + mag_err_list, + flux[mag_item], + error[mag_item]]) # (mag), (mag), (W m-2 um-1), (W m-2 um-1) - h5_file.create_dataset(f'objects/{object_name}/'+item, - data=data) + dset = h5_file.create_dataset(f'objects/{object_name}/{mag_item}', + data=data) + + dset.attrs['n_phot'] = n_phot if spectrum is not None: read_spec = {} @@ -1264,19 +1334,16 @@ def get_object(self, if filters: for item in filters: - data = dset[item] - - magnitude[item] = np.asarray(data[0:2]) - flux[item] = np.asarray(data[2:4]) + magnitude[item] = dset[item][0:2] + flux[item] = dset[item][2:4] else: for key in dset.keys(): if key not in ['distance', 'spectrum']: for item in dset[key]: - name = key+'/'+item - - magnitude[name] = np.asarray(dset[name][0:2]) - flux[name] = np.asarray(dset[name][2:4]) + name = f'{key}/{item}' + magnitude[name] = dset[name][0:2] + flux[name] = dset[name][2:4] filters = list(magnitude.keys()) diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index cd242486..2338ba5d 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -25,7 +25,7 @@ def plot_spectrum(boxes, filters=None, residuals=None, - colors=None, + plot_kwargs=None, xlim=None, ylim=None, ylim_res=None, @@ -46,8 +46,8 @@ def plot_spectrum(boxes, Filter IDs for which the transmission profile is plotted. Not plotted if set to None. residuals : species.core.box.ResidualsBox, None Box with residuals of a fit. Not plotted if set to None. - colors : list(str, ), None - Colors to be used for the different boxes. Note that a box with residuals requires a tuple + plot_kwargs : list(str, ), None + TODO Colors to be used for the different boxes. Note that a box with residuals requires a tuple with two colors (i.e., for the photometry and spectrum). Automatic colors are used if set to None. xlim : tuple(float, float) @@ -84,9 +84,12 @@ def plot_spectrum(boxes, marker = itertools.cycle(('o', 's', '*', 'p', '<', '>', 'P', 'v', '^')) - if colors is not None and len(boxes) != len(colors): + if plot_kwargs is None: + plot_kwargs = [] + + elif plot_kwargs is not None and len(boxes) != len(plot_kwargs): raise ValueError(f'The number of \'boxes\' ({len(boxes)}) should be equal to the ' - f'number of \'colors\' ({len(colors)}).') + f'number of items in \'plot_kwargs\' ({len(plot_kwargs)}).') if residuals is not None and filters is not None: plt.figure(1, figsize=figsize) @@ -253,10 +256,10 @@ def plot_spectrum(boxes, if residuals is not None: ax3.set_xscale(scale[0]) - color_obj_phot = None - color_obj_spec = None - for j, boxitem in enumerate(boxes): + if j < len(boxes): + plot_kwargs.append(None) + if isinstance(boxitem, (box.SpectrumBox, box.ModelBox)): wavelength = boxitem.wavelength flux = boxitem.flux @@ -317,9 +320,8 @@ def plot_spectrum(boxes, else: label = None - if colors: - ax1.plot(wavelength, masked/scaling, color=colors[j], lw=0.5, - label=label, zorder=2) + if plot_kwargs[j]: + ax1.plot(wavelength, masked/scaling, zorder=2, label=label, **plot_kwargs[j]) else: ax1.plot(wavelength, masked/scaling, lw=0.5, label=label, zorder=2) @@ -343,27 +345,23 @@ def plot_spectrum(boxes, data = np.array(flux, dtype=np.float64) masked = np.ma.array(data, mask=np.isnan(data)) - if colors: - ax1.plot(wavelength, masked/scaling, lw=0.2, color=colors[j], - alpha=0.5, zorder=1) + if plot_kwargs[j]: + ax1.plot(wavelength, masked/scaling, zorder=1, **plot_kwargs[j]) else: ax1.plot(wavelength, masked/scaling, lw=0.2, alpha=0.5, zorder=1) elif isinstance(boxitem, box.PhotometryBox): - marker = next(marker) - for i, item in enumerate(boxitem.wavelength): transmission = read_filter.ReadFilter(boxitem.filter_name[i]) fwhm = transmission.filter_fwhm() - if colors: + if plot_kwargs[j]: ax1.errorbar(item, boxitem.flux[i][0]/scaling, xerr=fwhm/2., - yerr=boxitem.flux[i][1], marker=marker, ms=6, color=colors[j], - zorder=3) + yerr=boxitem.flux[i][1]/scaling, zorder=3, **plot_kwargs[j]) else: ax1.errorbar(item, boxitem.flux[i][0]/scaling, xerr=fwhm/2., - yerr=boxitem.flux[i][1], marker=marker, ms=6, color='black', - zorder=3) + yerr=boxitem.flux[i][1]/scaling, marker='s', ms=6, color='black', + zorder=3) elif isinstance(boxitem, box.ObjectBox): if boxitem.flux is not None: @@ -372,48 +370,77 @@ def plot_spectrum(boxes, wavelength = transmission.mean_wavelength() fwhm = transmission.filter_fwhm() - if colors is None: - ax1.errorbar(wavelength, boxitem.flux[item][0]/scaling, xerr=fwhm/2., - yerr=boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=3, - markerfacecolor=color_obj_phot) + if not plot_kwargs[j] or item not in plot_kwargs[j]: + if not plot_kwargs[j]: + plot_kwargs[j] = {} + + if isinstance(boxitem.flux[item][0], np.ndarray): + for i in range(boxitem.flux[item].shape[1]): + + plot_obj = ax1.errorbar(wavelength, boxitem.flux[item][0, i]/scaling, xerr=fwhm/2., + yerr=boxitem.flux[item][1, i]/scaling, marker='s', ms=5, zorder=3) + + else: + plot_obj = ax1.errorbar(wavelength, boxitem.flux[item][0]/scaling, xerr=fwhm/2., + yerr=boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=3) + + plot_kwargs[j][item] = {'marker': 's', 'ms': 5., 'color': plot_obj[0].get_color()} else: - color_obj_phot = colors[j][0] + if isinstance(boxitem.flux[item][0], np.ndarray): + if not isinstance(plot_kwargs[j][item], list): + raise ValueError(f'A list with {boxitem.flux[item].shape[1]} ' + f'dictionaries are required because the filter ' + f'{item} has {boxitem.flux[item].shape[1]} ' + f'values.') + + for i in range(boxitem.flux[item].shape[1]): - ax1.errorbar(wavelength, boxitem.flux[item][0]/scaling, xerr=fwhm/2., - yerr=boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=3, - color=color_obj_phot, markerfacecolor=color_obj_phot) + ax1.errorbar(wavelength, boxitem.flux[item][0, i]/scaling, xerr=fwhm/2., + yerr=boxitem.flux[item][1, i]/scaling, zorder=3, **plot_kwargs[j][item][i]) + + else: + ax1.errorbar(wavelength, boxitem.flux[item][0]/scaling, xerr=fwhm/2., + yerr=boxitem.flux[item][1]/scaling, zorder=3, **plot_kwargs[j][item]) if boxitem.spectrum is not None: for key, value in boxitem.spectrum.items(): masked = np.ma.array(boxitem.spectrum[key][0], mask=np.isnan(boxitem.spectrum[key][0])) - if colors is None: + if not plot_kwargs[j] or key not in plot_kwargs[j]: ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling, ms=2, marker='s', zorder=2.5, ls='none') else: - color_obj_spec = colors[j][1] - ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling, - marker='o', ms=2, zorder=2.5, color=color_obj_spec, - markerfacecolor=color_obj_spec, ls='none') + zorder=2.5, **plot_kwargs[j][key]) elif isinstance(boxitem, box.SynphotBox): + for i, find_item in enumerate(boxes): + if isinstance(find_item, box.ObjectBox): + obj_index = i + break + for item in boxitem.flux: transmission = read_filter.ReadFilter(item) wavelength = transmission.mean_wavelength() fwhm = transmission.filter_fwhm() - if colors is None: + if not plot_kwargs[obj_index] or item not in plot_kwargs[obj_index]: ax1.errorbar(wavelength, boxitem.flux[item]/scaling, xerr=fwhm/2., yerr=None, - alpha=0.7, marker='s', ms=5, zorder=4, markerfacecolor='white') + alpha=0.7, marker='s', ms=5, zorder=4, mfc='white') else: - ax1.errorbar(wavelength, boxitem.flux[item]/scaling, xerr=fwhm/2., yerr=None, - alpha=0.7, marker='s', ms=5, zorder=4, color=colors[j], - markerfacecolor='white') + if isinstance(plot_kwargs[obj_index][item], list): + # In case of multiple photometry values for the same filter, use the + # plot_kwargs of the first data point + ax1.errorbar(wavelength, boxitem.flux[item]/scaling, xerr=fwhm/2., yerr=None, + zorder=4, mfc='white', **plot_kwargs[obj_index][item][0]) + + else: + ax1.errorbar(wavelength, boxitem.flux[item]/scaling, xerr=fwhm/2., yerr=None, + zorder=4, mfc='white', **plot_kwargs[obj_index][item]) if filters is not None: for i, item in enumerate(filters): @@ -423,23 +450,38 @@ def plot_spectrum(boxes, ax2.plot(data[0, ], data[1, ], '-', lw=0.7, color='black', zorder=1) if residuals is not None: + for i, find_item in enumerate(boxes): + if isinstance(find_item, box.ObjectBox): + obj_index = i + break + res_max = 0. if residuals.photometry is not None: - ax3.plot(residuals.photometry[0, ], residuals.photometry[1, ], marker='s', - ms=5, linestyle='none', color=color_obj_phot, zorder=2) + for item in residuals.photometry: + if not plot_kwargs[obj_index] or item not in plot_kwargs[obj_index]: + ax3.plot(residuals.photometry[item][0], residuals.photometry[item][1], marker='s', + ms=5, linestyle='none', zorder=2) + + else: + if residuals.photometry[item].ndim == 1: + ax3.plot(residuals.photometry[item][0], residuals.photometry[item][1], zorder=2, + **plot_kwargs[obj_index][item]) + + elif residuals.photometry[item].ndim == 2: + for i in range(residuals.photometry[item].shape[1]): + ax3.plot(residuals.photometry[item][0, i], residuals.photometry[item][1, i], zorder=2, + **plot_kwargs[obj_index][item][i]) - res_max = np.nanmax(np.abs(residuals.photometry[1, ])) + res_max = np.nanmax(np.abs(residuals.photometry[item][1])) if residuals.spectrum is not None: for key, value in residuals.spectrum.items(): - if colors is None: - ax3.plot(value[:, 0], value[:, 1], marker='o', ms=2, - linestyle='none', zorder=1) + if not plot_kwargs[obj_index] or key not in plot_kwargs[obj_index]: + ax3.plot(value[:, 0], value[:, 1], marker='o', ms=2, ls='none', zorder=1) else: - ax3.plot(value[:, 0], value[:, 1], marker='o', - ms=2, linestyle='none', color=color_obj_spec, zorder=1) + ax3.plot(value[:, 0], value[:, 1], zorder=1, **plot_kwargs[obj_index][key]) max_tmp = np.nanmax(np.abs(value[:, 1])) diff --git a/species/util/phot_util.py b/species/util/phot_util.py index 9c1ca700..e8adf8b3 100644 --- a/species/util/phot_util.py +++ b/species/util/phot_util.py @@ -5,8 +5,6 @@ import math import warnings -import matplotlib.pyplot as plt - import spectres import numpy as np @@ -180,13 +178,22 @@ def get_residuals(datatype, filters=filters, parameters=parameters) - res_phot = np.zeros((2, len(objectbox.flux))) + res_phot = {} - for i, item in enumerate(filters): + for item in filters: transmission = read_filter.ReadFilter(item) + res_phot[item] = np.zeros(objectbox.flux[item].shape) + + if objectbox.flux[item].ndim == 1: + res_phot[item][0] = transmission.mean_wavelength() + res_phot[item][1] = (objectbox.flux[item][0]-model_phot.flux[item]) / \ + objectbox.flux[item][1] - res_phot[0, i] = transmission.mean_wavelength() - res_phot[1, i] = (objectbox.flux[item][0]-model_phot.flux[item])/objectbox.flux[item][1] + elif objectbox.flux[item].ndim == 2: + for j in range(objectbox.flux[item].shape[1]): + res_phot[item][0, j] = transmission.mean_wavelength() + res_phot[item][1, j] = (objectbox.flux[item][0, j]-model_phot.flux[item]) / \ + objectbox.flux[item][1, j] else: res_phot = None @@ -250,8 +257,13 @@ def get_residuals(datatype, print('Residuals (sigma):') if res_phot is not None: - for i, item in enumerate(filters): - print(f' - {item}: {res_phot[1, i]:.2f}') + for item in filters: + if res_phot[item].ndim == 1: + print(f' - {item}: {res_phot[item][1]:.2f}') + + elif res_phot[item].ndim == 2: + for j in range(res_phot[item].shape[1]): + print(f' - {item}: {res_phot[item][1, j]:.2f}') if res_spec is not None: for key in objectbox.spectrum: diff --git a/test/test_read/test_calibration.py b/test/test_read/test_calibration.py index ef5017cb..5f550a71 100644 --- a/test/test_read/test_calibration.py +++ b/test/test_read/test_calibration.py @@ -42,8 +42,8 @@ def test_get_spectrum(self): read_calib = species.ReadCalibration('vega', filter_name='Paranal/NACO.Lp') spec_box = read_calib.get_spectrum(self.model_param, apply_mask=True, spec_res=100.) - assert np.sum(spec_box.wavelength) == pytest.approx(175.34279961519331, rel=self.limit, abs=0.) - assert np.sum(spec_box.flux) == pytest.approx(2.3081127325020912e-09, rel=self.limit, abs=0.) + assert np.sum(spec_box.wavelength) == pytest.approx(183.34793783628004, rel=self.limit, abs=0.) + assert np.sum(spec_box.flux) == pytest.approx(2.363239096767195e-09, rel=self.limit, abs=0.) def test_get_flux(self): read_calib = species.ReadCalibration('vega', filter_name='Paranal/NACO.H') diff --git a/test/test_read/test_isochrone.py b/test/test_read/test_isochrone.py index 3a12a688..805fc696 100644 --- a/test/test_read/test_isochrone.py +++ b/test/test_read/test_isochrone.py @@ -80,10 +80,10 @@ def test_get_color_magnitude(self): assert colormag_box.color.shape == (10, ) assert colormag_box.magnitude.shape == (10, ) - assert np.sum(colormag_box.color) == pytest.approx(2.483071377517346, + assert np.sum(colormag_box.color) == pytest.approx(2.496406543555395, rel=self.limit, abs=0.) - assert np.sum(colormag_box.magnitude) == pytest.approx(109.57418536215742, + assert np.sum(colormag_box.magnitude) == pytest.approx(109.59296195800431, rel=self.limit, abs=0.) assert np.sum(colormag_box.sptype) == pytest.approx(400., rel=self.limit, abs=0.) @@ -101,10 +101,10 @@ def test_get_color_color(self): assert colorcolor_box.color1.shape == (10, ) assert colorcolor_box.color2.shape == (10, ) - assert np.sum(colorcolor_box.color1) == pytest.approx(2.483071377517346, + assert np.sum(colorcolor_box.color1) == pytest.approx(2.496406543555395, rel=self.limit, abs=0.) - assert np.sum(colorcolor_box.color2) == pytest.approx(3.3474795384424585, + assert np.sum(colorcolor_box.color2) == pytest.approx(3.3537987656436457, rel=self.limit, abs=0.) assert np.sum(colorcolor_box.sptype) == pytest.approx(400., rel=self.limit, abs=0.) diff --git a/test/test_read/test_model.py b/test/test_read/test_model.py index 1fee0bfe..023e130b 100644 --- a/test/test_read/test_model.py +++ b/test/test_read/test_model.py @@ -43,36 +43,36 @@ def test_get_model(self): magnitude=False, smooth=True) - assert np.sum(model_box.wavelength) == pytest.approx(82.34125130167924, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(1.5427506511556731e-12, rel=self.limit, abs=0.) + assert np.sum(model_box.wavelength) == pytest.approx(83.77689756715003, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(1.5697815415870652e-12, rel=self.limit, abs=0.) model_box = read_model.get_model(self.model_param, spec_res=100., magnitude=True, smooth=True) - assert np.sum(model_box.wavelength) == pytest.approx(82.34125130167924, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(576.6670455683387, rel=self.limit, abs=0.) + assert np.sum(model_box.wavelength) == pytest.approx(83.77689756715003, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(588.6995381002678, rel=self.limit, abs=0.) def test_get_data(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') model_box = read_model.get_data(self.model_param) - assert np.sum(model_box.wavelength) == pytest.approx(90.85089938365051, rel=self.limit, abs=0.) - assert np.sum(model_box.flux) == pytest.approx(1.6100966071681859e-12, rel=self.limit, abs=0.) + assert np.sum(model_box.wavelength) == pytest.approx(92.26773310928259, rel=self.limit, abs=0.) + assert np.sum(model_box.flux) == pytest.approx(1.6346965666899605e-12, rel=self.limit, abs=0.) def test_get_flux(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') flux = read_model.get_flux(self.model_param) - assert flux[0] == pytest.approx(3.3369304657212616e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(3.336845306372948e-14, rel=self.limit, abs=0.) def test_get_magnitude(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') magnitude = read_model.get_magnitude(self.model_param) - assert magnitude[0] == pytest.approx(11.357113310356498, rel=self.limit, abs=0.) - assert magnitude[1] == pytest.approx(11.357113310356498, rel=self.limit, abs=0.) + assert magnitude[0] == pytest.approx(11.357141018985255, rel=self.limit, abs=0.) + assert magnitude[1] == pytest.approx(11.357141018985255, rel=self.limit, abs=0.) def test_get_bounds(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') @@ -85,7 +85,7 @@ def test_get_wavelengths(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') wavelengths = read_model.get_wavelengths() - assert np.sum(wavelengths) == pytest.approx(801.5391789533717, rel=1e-7, abs=0.) + assert np.sum(wavelengths) == pytest.approx(813.2224003071026, rel=1e-7, abs=0.) def test_get_points(self): read_model = species.ReadModel('ames-cond', filter_name='Paranal/NACO.H') diff --git a/test/test_read/test_planck.py b/test/test_read/test_planck.py index 2f5f695c..49190d8a 100644 --- a/test/test_read/test_planck.py +++ b/test/test_read/test_planck.py @@ -38,15 +38,15 @@ def test_get_spectrum(self): assert modelbox.wavelength.shape == (42, ) assert modelbox.flux.shape == (42, ) - assert np.sum(modelbox.wavelength) == pytest.approx(52.7026751397061, rel=self.limit, abs=0.) - assert np.sum(modelbox.flux) == pytest.approx(8.332973825446955e-13, rel=self.limit, abs=0.) + assert np.sum(modelbox.wavelength) == pytest.approx(52.581085755096346, rel=self.limit, abs=0.) + assert np.sum(modelbox.flux) == pytest.approx(8.322089100757744e-13, rel=self.limit, abs=0.) def test_get_flux(self): read_planck = species.ReadPlanck(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}) - assert flux[0] == pytest.approx(1.9888949873704357e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(1.9888812770713738e-14, rel=self.limit, abs=0.) synphot = species.SyntheticPhotometry(filter_name='MKO/NSFCam.J') flux = read_planck.get_flux({'teff': 2000., 'radius': 1., 'distance': 10.}, synphot=synphot) - assert flux[0] == pytest.approx(1.9888949873704357e-14, rel=self.limit, abs=0.) + assert flux[0] == pytest.approx(1.9888812770713738e-14, rel=self.limit, abs=0.)