Skip to content

Commit

Permalink
Object support for duplicate filter names
Browse files Browse the repository at this point in the history
  • Loading branch information
Tomas Stolker committed Apr 30, 2020
1 parent 1eec252 commit dd7566d
Show file tree
Hide file tree
Showing 10 changed files with 270 additions and 159 deletions.
47 changes: 16 additions & 31 deletions species/analysis/photometry.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]

Expand All @@ -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.')

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
21 changes: 14 additions & 7 deletions species/core/box.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
6 changes: 2 additions & 4 deletions species/data/companions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
151 changes: 109 additions & 42 deletions species/data/database.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
-------
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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 = {}
Expand Down Expand Up @@ -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())

Expand Down
Loading

0 comments on commit dd7566d

Please sign in to comment.