From 1041eb5fbb52acbdfecea7276baf37f5ffadca96 Mon Sep 17 00:00:00 2001 From: Tomas Stolker Date: Tue, 19 Mar 2019 17:13:32 +0100 Subject: [PATCH] improved interpolation --- species/__init__.py | 9 +- species/analysis/fit_model.py | 285 +++++++++++++++--------- species/analysis/photometry.py | 176 +++++++-------- species/core/box.py | 68 +++++- species/data/companions.py | 8 +- species/data/database.py | 383 +++++++++++++++++++++------------ species/plot/plot_mcmc.py | 40 +++- species/plot/plot_spectrum.py | 112 +++++++--- species/read/read_model.py | 382 ++++++++++++++++---------------- species/read/read_object.py | 81 +++++-- species/util/data_util.py | 83 ++++--- species/util/phot_util.py | 133 ++++++++++++ species/util/plot_util.py | 2 +- species/util/read_util.py | 2 +- 14 files changed, 1133 insertions(+), 631 deletions(-) create mode 100644 species/util/phot_util.py diff --git a/species/__init__.py b/species/__init__.py index 2de9b62d..da398d1c 100644 --- a/species/__init__.py +++ b/species/__init__.py @@ -2,10 +2,7 @@ from species.analysis.fit_spectrum import FitSpectrum -from species.analysis.photometry import SyntheticPhotometry, \ - apparent_to_absolute, \ - multi_photometry - +from species.analysis.photometry import SyntheticPhotometry from species.read.read_calibration import ReadCalibration @@ -42,6 +39,10 @@ from species.plot.plot_spectrum import plot_spectrum +from species.util.phot_util import apparent_to_absolute, \ + multi_photometry, \ + get_residuals + from species.util.read_util import get_mass, \ add_luminosity diff --git a/species/analysis/fit_model.py b/species/analysis/fit_model.py index 102ba270..36e497f5 100644 --- a/species/analysis/fit_model.py +++ b/species/analysis/fit_model.py @@ -1,17 +1,19 @@ """ -Text +Module for fitting atmospheric models. """ import sys import math import emcee +import spectres import progress.bar import numpy as np from species.analysis import photometry from species.data import database from species.read import read_model, read_object +from species.util import read_util MIN_CHISQ = np.inf @@ -23,19 +25,25 @@ def lnprior(param, modelpar, prior): """ - :param param: Parameter values. - :type param: numpy.ndarray - :param bounds: Parameter boundaries. - :type bounds: dict - :param modelpar: Parameter names. - :type modelpar: tuple(str, ) - :param prior: Gaussian prior on one of the parameters. Currently only possible for the mass, - e.g. ('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of - 3 Mjup. Not used if set to None. - :type prior: tuple(str, float, float) - - :return: Log prior probability. - :rtype: float + Function for the prior probability. + + Parameters + ---------- + param : numpy.ndarray + Parameter values. + bounds : dict + Parameter boundaries. + modelpar : list(str, ) + Parameter names. + prior : tuple(str, float, float) + Gaussian prior on one of the parameters. Currently only possible for the mass, e.g. + ('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. Not + used if set to None. + + Returns + ------- + float + Log prior probability. """ if prior: @@ -52,7 +60,7 @@ def lnprior(param, ln_prior = 0. elif prior[0] == 'mass': - mass = read_model.get_mass(modeldict) + mass = read_util.get_mass(modeldict) ln_prior = -0.5*(mass-prior[1])**2/prior[2]**2 else: @@ -67,26 +75,34 @@ def lnlike(param, modelphot, objphot, synphot, - sampling, - distance): + distance, + spectrum, + instrument, + modelspec): """ - :param param: - :type param: - :param modelpar: - :type modelpar: - :param modelphot: - :type modelphot: - :param objphot: - :type objphot: - :param synphot: - :type synphot: - :param sampling: - :type sampling: - :param distance: - :type distance: - - :return: Log likelihood probability. - :rtype: float + Function for the likelihood probability. + + Parameters + ---------- + param : numpy.ndarray + Parameter values. + modelpar : list(str, ) + Parameter names. + modelphot : list('species.read.read_model.ReadModel, ) + objphot : list(tuple(float, float), ) + synphot : list(species.analysis.photometry.SyntheticPhotometry, ) + distance : float + Distance (pc). + spectrum : numpy.ndarray + Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). + instrument : str + Instrument that was used for the spectrum (currently only 'gpi' possible). + modelspec : species.read.read_model.ReadModel + + Returns + ------- + float + Log likelihood probability. """ global MIN_CHISQ @@ -100,9 +116,20 @@ def lnlike(param, chisq = 0. - for i, item in enumerate(objphot): - flux = modelphot[i].get_photometry(paramdict, sampling, synphot[i]) - chisq += (item[0]-flux)**2 / item[1]**2 + if objphot is not None: + for i, item in enumerate(objphot): + flux = modelphot[i].get_photometry(paramdict, synphot[i]) + chisq += (item[0]-flux)**2 / item[1]**2 + + if spectrum is not None: + model = modelspec.get_model(paramdict, ('original', (0.9, 2.5))) + + flux_new = spectres.spectres(new_spec_wavs=spectrum[:, 0], + old_spec_wavs=model.wavelength, + spec_fluxes=model.flux, + spec_errs=None) + + chisq += np.nansum((spectrum[:, 1] - flux_new)**2/spectrum[:, 2]**2) if chisq < MIN_CHISQ: MIN_CHISQ = chisq @@ -117,31 +144,41 @@ def lnprob(param, modelphot, objphot, synphot, - sampling, distance, - prior): + prior, + spectrum, + instrument, + modelspec): """ - :param param: - :type param: - :param bounds: - :type bounds: - :param modelpar: - :type modelpar: - :param modelphot: - :type modelphot: - :param objphot: - :type objphot: - :param synphot: - :type synphot: - :param sampling: - :type sampling: - :param distance: - :type distance: - :param prior: Gaussian prior. Not used if set to None. - :type prior: tuple(str, float, float) - - :return: - :rtype: + Function for the posterior probability. + + Parameters + ---------- + param : numpy.ndarray + Parameter values. + bounds : dict + Parameter boundaries. + modelpar : list(str, ) + Parameter names. + modelphot : list('species.read.read_model.ReadModel, ) + objphot : list(tuple(float, float), ) + synphot : list(species.analysis.photometry.SyntheticPhotometry, ) + distance : float + Distance (pc). + prior : tuple(str, float, float) + Gaussian prior on one of the parameters. Currently only possible for the mass, e.g. + ('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. Not + used if set to None. + spectrum : numpy.ndarray + Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). + instrument : str + Instrument that was used for the spectrum (currently only 'gpi' possible). + modelspec : species.read.read_model.ReadModel + + Returns + ------- + float + Log posterior probability. """ ln_prior = lnprior(param, bounds, modelpar, prior) @@ -155,53 +192,57 @@ def lnprob(param, modelphot, objphot, synphot, - sampling, - distance) + distance, + spectrum, + instrument, + modelspec) return ln_prob class FitModel: """ - Text + Fit atmospheric model spectra to photometric data. """ def __init__(self, objname, filters, model, - sampling, - bounds): + bounds, + data): """ - :param objname: Object name in the database. - :type objname: str - :param filters: Filter IDs for which the photometry is selected. All available - photometry of the object is selected if set to None. - :type filters: tuple(str, ) - :name model: Atmospheric model. - :type model: str - :name sampling: Wavelength sampling for the computation of synthetic photometry - ('specres' or 'gaussian'). - :type sampling: tuple - :name bounds: Parameter boundaries. Full parameter range is used if None or not specified. - The radius parameter range is set to 0-5 Rjup if not specified. - :type bounds: dict - - :return: None + Parameters + ---------- + objname : str + Object name in the database. + filters : tuple(str, ) + Filter IDs for which the photometry is selected. All available photometry of the + object is selected if set to None. + model : str + Atmospheric model. + bounds : dict + Parameter boundaries. Full parameter range is used if None or not specified. The + radius parameter range is set to 0-5 Rjup if not specified. + data : tuple(bool, bool) + Data to use for the fit ('photometry' and/or 'spectrum'). The first boolean sets + the inclusion of photometry and the second boolean the inclusion of a spectrum. + + Returns + ------- + None """ self.object = read_object.ReadObject(objname) self.distance = self.object.get_distance() self.model = model - self.sampling = sampling self.bounds = bounds - self.objphot = [] - self.modelphot = [] - self.synphot = [] + if not data[0] and not data[1]: + raise ValueError('No photometric or spectral data has been selected.') - if self.bounds and 'teff' in self.bounds: + if self.bounds is not None and 'teff' in self.bounds: teff_bound = self.bounds['teff'] else: teff_bound = None @@ -221,21 +262,40 @@ def __init__(self, if 'radius' not in self.bounds: self.bounds['radius'] = (0., 5.) - if filters is None: - species_db = database.Database() - objectbox = species_db.get_object(objname, None) - filters = objectbox.filter + if data[0]: + self.objphot = [] + self.modelphot = [] + self.synphot = [] - for item in filters: - readmodel = read_model.ReadModel(self.model, item, teff_bound) - readmodel.interpolate() - self.modelphot.append(readmodel) + if not filters: + species_db = database.Database() + objectbox = species_db.get_object(objname, None) + filters = objectbox.filter - sphot = photometry.SyntheticPhotometry(item) - self.synphot.append(sphot) + for item in filters: + readmodel = read_model.ReadModel(self.model, item, teff_bound) + readmodel.interpolate() + self.modelphot.append(readmodel) - obj_phot = self.object.get_photometry(item) - self.objphot.append((obj_phot[2], obj_phot[3])) + sphot = photometry.SyntheticPhotometry(item) + self.synphot.append(sphot) + + obj_phot = self.object.get_photometry(item) + self.objphot.append((obj_phot[2], obj_phot[3])) + + else: + self.objphot = None + self.modelphot = None + self.synphot = None + + if data[1]: + self.spectrum = self.object.get_spectrum() + self.instrument = self.object.get_instrument() + self.modelspec = read_model.ReadModel(self.model, (0.9, 2.5), teff_bound) + else: + self.spectrum = None + self.instrument = None + self.modelspec = None self.modelpar = readmodel.get_parameters() self.modelpar.append('radius') @@ -248,7 +308,30 @@ def run_mcmc(self, prior=None, ncpu=1): """ - :return: None + Function to run the MCMC sampler. + + Parameters + ---------- + nwalkers : int + Number of walkers. + nsteps : int + Number of steps per walker. + guess : dict + Guess for the parameter values. Random values between the boundary values are used + if set to None. + tag : str + Database tag where the MCMC samples are stored. + prior : tuple(str, float, float) + Gaussian prior on one of the parameters. Currently only possible for the mass, e.g. + ('mass', 13., 3.) for an expected mass of 13 Mjup with an uncertainty of 3 Mjup. Not + used if set to None. + ncpu : int + Number of parallel processes. Due to the additional overhead, a value larger than 1 + will not necessarily result in a decrease in computation time. + + Returns + ------- + None """ global MIN_CHISQ @@ -280,9 +363,11 @@ def run_mcmc(self, self.modelphot, self.objphot, self.synphot, - self.sampling, self.distance, - prior]), + prior, + self.spectrum, + self.instrument, + self.modelspec]), threads=ncpu) progbar = progress.bar.Bar('\rRunning MCMC...', diff --git a/species/analysis/photometry.py b/species/analysis/photometry.py index b706a4f4..4b412eb6 100644 --- a/species/analysis/photometry.py +++ b/species/analysis/photometry.py @@ -1,5 +1,5 @@ """ -Photometry module. +Module for creating synthetic photometry. """ import os @@ -9,70 +9,25 @@ import numpy as np from species.data import database -from species.core import box -from species.read import read_filter, read_calibration, read_model - - -def multi_photometry(datatype, - spectrum, - filters, - parameters): - """ - :param datatype: Data type ('model' or 'calibration'). - :type datatype: str - :param spectrum: Spectrum name (e.g., 'drift-phoenix'). - :type spectrum: str - :param filters: Filter IDs. - :type filters: tuple(str, ) - :param parameters: Model parameter values. - :type parameters: dict - - :return: Box with synthetic photometry. - :rtype: species.core.box.SynphotBox - """ - - flux = {} - - if datatype == 'model': - for item in filters: - readmodel = read_model.ReadModel(spectrum, item) - flux[item] = readmodel.get_photometry(parameters, ('specres', 100.)) - - elif datatype == 'calibration': - for item in filters: - readcalib = read_calibration.ReadCalibration(spectrum, item) - flux[item] = readcalib.get_photometry(parameters) - - return box.create_box('synphot', name='synphot', flux=flux) - - -def apparent_to_absolute(app_mag, - distance): - """ - :param app_mag: Apparent magnitude (mag). - :type app_mag: float or numpy.ndarray - :param distance: Distance (pc). - :type distance: float or numpy.ndarray - - :return: Absolute magnitude (mag). - :rtype: float or numpy.ndarray - """ - - return app_mag - 5.*np.log10(distance) + 5. +from species.read import read_filter, read_calibration class SyntheticPhotometry: """ - Text + Create synthetic photometry from a spectrum. """ def __init__(self, filter_name): """ - :param filter_name: Filter name. - :type filter_name: str - - :return: None + Parameters + ---------- + filter_name : str + Filter ID. + + Returns + ------- + None """ self.filter_name = filter_name @@ -89,12 +44,15 @@ def __init__(self, def zero_point(self, wl_range): """ - :param wl_range: Wavelength range (micron). The range from the filter transmission - curve is used if set to None. - :type wl_range: float - - :return: tuple(float, float) - :rtype: + Parameters + ---------- + wl_range : float + Wavelength range (micron). The range from the filter transmission curve is used if + set to None. + + Returns + ------- + tuple(float, float) """ if wl_range is None: @@ -129,13 +87,17 @@ def spectrum_to_photometry(self, wavelength, flux_density): """ - :param wavelength: Wavelength (micron). - :type wavelength: numpy.ndarray - :param flux_density: Flux density (W m-2 micron-1). - :type flux_density: numpy.ndarray - - :return: Average flux density (W m-2 micron-1). - :rtype: float or numpy.ndarray + Parameters + ---------- + wavelength : numpy.ndarray + Wavelength (micron). + flux_density : numpy.ndarray + Flux density (W m-2 micron-1). + + Returns + ------- + float or numpy.ndarray + Average flux density (W m-2 micron-1). """ if not self.filter_interp: @@ -147,6 +109,10 @@ def spectrum_to_photometry(self, indices = np.where((self.wl_range[0] < wavelength) & (wavelength < self.wl_range[1]))[0] + if indices.size == 1: + raise ValueError("Calculating synthetic photometry requires more than one " + "wavelength point.") + wavelength = wavelength[indices] flux_density = flux_density[indices] @@ -169,6 +135,10 @@ def spectrum_to_photometry(self, indices = np.where((self.wl_range[0] <= wavelength) & (wavelength <= self.wl_range[1]))[0] + if indices.size == 1: + raise ValueError("Calculating synthetic photometry requires more than one " + "wavelength point.") + wavelength = wavelength[indices] flux_density = flux_density[indices] @@ -194,15 +164,21 @@ def spectrum_to_magnitude(self, flux_density, distance=None): """ - :param wavelength: Wavelength (micron). - :type wavelength: numpy.ndarray - :param flux_density: Flux density (W m-2 micron-1). - :type flux_density: numpy.ndarray - :param distance: Distance (pc). No absolute magnitude is calculated if set to None. - :type distance: float - - :return: Flux (W m-2). - :rtype: float or numpy.ndarray + Parameters + ---------- + wavelength : numpy.ndarray + Wavelength (micron). + flux_density : numpy.ndarray + Flux density (W m-2 micron-1). + distance : float + Distance (pc). No absolute magnitude is calculated if set to None. + + Returns + ------- + float or numpy.ndarray + Apparent magnitude (mag). + float or numpy.ndarray + Absolute magnitude (mag). """ vega_mag = 0.03 # [mag] @@ -231,15 +207,21 @@ def magnitude_to_flux(self, error, zp_flux=None): """ - :param magnitude: Magnitude (mag). - :type magnitude: float - :param error: Error (mag). Not used if set to None. - :type error: float - :param zp_flux: Zero point flux (W m-2 micron-1). Calculated if set to None. - :type zp_flux: float - - :return: Flux (W m-2 micron-1), lower error, upper error - :rtype: float, tuple(float, float) + Parameters + ---------- + magnitude : float + Magnitude (mag). + error : float + Error (mag). Not used if set to None. + zp_flux : float + Zero point flux (W m-2 micron-1). Calculated if set to None. + + Returns + ------- + float + Flux density (W m-2 micron-1). + float + Error (W m-2 micron-1). """ vega_mag = 0.03 # [mag] @@ -263,13 +245,19 @@ def flux_to_magnitude(self, flux, distance): """ - :param flux: Flux density (W m-2 micron-1). - :type flux: float - :param error: Distance (pc). - :type error: float - - :return: Apparent magnitude (mag), absolute magnitude (mag). - :rtype: float, float + Parameters + ---------- + flux : float + Flux density (W m-2 micron-1). + error : float + Distance (pc). + + Returns + ------- + float + Apparent magnitude (mag). + float + Absolute magnitude (mag). """ vega_mag = 0.03 # [mag] diff --git a/species/core/box.py b/species/core/box.py index ab26cff1..db0cc6a8 100644 --- a/species/core/box.py +++ b/species/core/box.py @@ -2,10 +2,14 @@ Box module. """ +import sys + def open_box(box): """ - return: + Returns + ------- + None """ for item in box.__dict__.keys(): @@ -14,7 +18,9 @@ def open_box(box): def create_box(boxtype, **kwargs): """ - :return: + Returns + ------- + species.core.box """ if boxtype == 'colormag': @@ -50,6 +56,7 @@ def create_box(boxtype, **kwargs): box.magnitude = kwargs['magnitude'] box.flux = kwargs['flux'] box.distance = kwargs['distance'] + box.spectrum = kwargs['spectrum'] elif boxtype == 'photometry': box = PhotometryBox() @@ -57,6 +64,12 @@ def create_box(boxtype, **kwargs): box.wavelength = kwargs['wavelength'] box.flux = kwargs['flux'] + elif boxtype == 'residuals': + box = ResidualsBox() + box.name = kwargs['name'] + box.photometry = kwargs['photometry'] + box.spectrum = kwargs['spectrum'] + elif boxtype == 'samples': box = SamplesBox() box.spectrum = kwargs['spectrum'] @@ -90,7 +103,9 @@ class ColorMagBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.library = None @@ -109,7 +124,9 @@ class ColorColorBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.library = None @@ -127,7 +144,9 @@ class ModelBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.model = None @@ -144,7 +163,9 @@ class ObjectBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.name = None @@ -152,6 +173,7 @@ def __init__(self): self.magnitude = None self.flux = None self.distance = None + self.spectrum = None class PhotometryBox: @@ -161,7 +183,9 @@ class PhotometryBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.name = None @@ -169,6 +193,23 @@ def __init__(self): self.flux = None +class ResidualsBox: + """ + Text + """ + + def __init__(self): + """ + Returns + ------- + None + """ + + self.name = None + self.photometry = None + self.spectrum = None + + class SamplesBox: """ Text @@ -176,7 +217,9 @@ class SamplesBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.spectrum = None @@ -192,12 +235,15 @@ class SpectrumBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.spectrum = None self.wavelength = None self.flux = None + self.error = None self.name = None self.simbad = None self.sptype = None @@ -211,7 +257,9 @@ class SynphotBox: def __init__(self): """ - :return: + Returns + ------- + None """ self.name = None diff --git a/species/data/companions.py b/species/data/companions.py index f3875c8e..041a26b4 100644 --- a/species/data/companions.py +++ b/species/data/companions.py @@ -1,11 +1,13 @@ """ -Text +Module with a function for data of directly imaged companions. """ def get_data(): """ - :return: Dictionary with the distances and apparent magnitudes of directly imaged companions. - :rtype: dict + Returns + ------- + dict + Dictionary with the distances and apparent magnitudes of directly imaged companions. """ data = {'beta Pic b':{'distance':19.75, diff --git a/species/data/database.py b/species/data/database.py index 8e193ca8..b09ac9d7 100644 --- a/species/data/database.py +++ b/species/data/database.py @@ -9,6 +9,7 @@ import h5py import emcee +import progress.bar import numpy as np from astropy.io import votable @@ -31,7 +32,9 @@ class Database: def __init__(self): """ - :return: None + Returns + ------- + None """ config_file = os.path.join(os.getcwd(), 'species_config.ini') @@ -44,7 +47,9 @@ def __init__(self): def list_items(self): """ - return: None + Returns + ------- + None """ sys.stdout.write('Database content:\n') @@ -52,12 +57,14 @@ def list_items(self): def descend(h5_object, seperator=''): """ - :param h5_object: - :type h5_object: h5py._hl.files.File, h5py._hl.group.Group, h5py._hl.dataset.Dataset - :param separator: - :type separator: str - - :return: None + Parameters + ---------- + h5_object : h5py._hl.files.File, h5py._hl.group.Group, h5py._hl.dataset.Dataset + separator : str + + Returns + ------- + None """ if isinstance(h5_object, (h5py._hl.files.File, h5py._hl.group.Group)): @@ -77,7 +84,9 @@ def descend(h5_object, def list_companions(self): """ - :return: None + Returns + ------- + None """ comp_phot = companions.get_data() @@ -90,10 +99,14 @@ def list_companions(self): def add_companion(self, name=None): """ - :param name: Companion name. All companions are added if set to None. - :type name: tuple(str, ) + Parameters + ---------- + tuple(str, ) + Companion name. All companions are added if set to None. - :return: None + Returns + ------- + None """ if isinstance(name, str): @@ -113,15 +126,18 @@ def add_filter(self, filter_id, filename=None): """ - :param filter_id: Filter ID from the SVO Filter Profile Service (e.g., 'Paranal/NACO.Lp'). - :type filter_id: str - :param filename: Filename with the filter profile. The first column should contain the - wavelength (micron) and the second column the transmission (no units). - The profile is downloaded from the SVO Filter Profile Service if set to - None. - :type filename: str + Parameters + ---------- + filter_id : str + Filter ID from the SVO Filter Profile Service (e.g., 'Paranal/NACO.Lp'). + filename : str + Filename with the filter profile. The first column should contain the wavelength + (micron) and the second column the transmission (no units). The profile is downloaded + from the SVO Filter Profile Service if set to None. - :return: None + Returns + ------- + None """ filter_split = filter_id.split('/') @@ -162,14 +178,18 @@ def add_model(self, wavelength=None, teff=None): """ - :param model: Model name. - :type model: str - :param wavelength: Wavelength (micron) range. - :type wavelength: tuple(float, float) - :param teff: Effective temperature (K) range. - :type teff: tuple(float, float) + Parameters + ---------- + model : str + Model name. + wavelength : tuple(float, float) + Wavelength (micron) range. + teff : tuple(float, float) + Effective temperature (K) range. - :return: None + Returns + ------- + None """ h5_file = h5py.File(self.database, 'a') @@ -189,29 +209,32 @@ def add_model(self, def add_object(self, object_name, - distance, - app_mag): - """ - :param object_name: Object name. - :type object_name: str - :param distance: Distance (pc). - :type distance: float - :param app_mag: Apparent magnitudes. - :type app_mag: dict - - :return: None + distance=None, + app_mag=None, + spectrum=None, + instrument=None): + """ + Parameters + ---------- + object_name: str + Object name. + distance : float + Distance (pc). Not written if set to None. + app_mag : dict + Apparent magnitudes. Not written if set to None. + spectrum : str + Spectrum filename. The first three columns should contain the wavelength (micron), + flux density (W m-2 micron-1), and the error (W m-2 micron-1). Not written if set + to None. + instrument : str + Instrument that was used for the spectrum (currently only 'gpi' possible). Not + used if set to None. + + Returns + ------- + None """ - flux = {} - error = {} - - for item in app_mag: - synphot = photometry.SyntheticPhotometry(item) - flux[item], error[item] = synphot.magnitude_to_flux(app_mag[item][0], app_mag[item][1]) - - sys.stdout.write('Adding object: '+object_name+'...') - sys.stdout.flush() - h5_file = h5py.File(self.database, 'a') if 'objects' not in h5_file: @@ -220,26 +243,54 @@ def add_object(self, if 'objects/'+object_name not in h5_file: h5_file.create_group('objects/'+object_name) - if 'objects/'+object_name+'/distance' in h5_file: - del h5_file['objects/'+object_name+'/distance'] + if distance: - h5_file.create_dataset('objects/'+object_name+'/distance', - data=distance, - dtype='f') # [pc] + if 'objects/'+object_name+'/distance' in h5_file: + del h5_file['objects/'+object_name+'/distance'] - for item in app_mag: - if 'objects/'+object_name+'/'+item in h5_file: - del h5_file['objects/'+object_name+'/'+item] + h5_file.create_dataset('objects/'+object_name+'/distance', + data=distance, + dtype='f') # [pc] - data = np.asarray([app_mag[item][0], - app_mag[item][1], - flux[item], - error[item]]) + if app_mag: - # [mag], [mag], [W m-2 micron-1], [W m-2 micron-1] - h5_file.create_dataset('objects/'+object_name+'/'+item, - data=data, - dtype='f') + flux = {} + error = {} + + for item in app_mag: + synphot = photometry.SyntheticPhotometry(item) + flux[item], error[item] = synphot.magnitude_to_flux(app_mag[item][0], + app_mag[item][1]) + + for item in app_mag: + if 'objects/'+object_name+'/'+item in h5_file: + del h5_file['objects/'+object_name+'/'+item] + + data = np.asarray([app_mag[item][0], + app_mag[item][1], + flux[item], + error[item]]) + + # [mag], [mag], [W m-2 micron-1], [W m-2 micron-1] + h5_file.create_dataset('objects/'+object_name+'/'+item, + data=data, + dtype='f') + + sys.stdout.write('Adding object: '+object_name+'...') + sys.stdout.flush() + + if spectrum: + + if 'objects/'+object_name+'/spectrum' in h5_file: + del h5_file['objects/'+object_name+'/spectrum'] + + data = np.loadtxt(spectrum) + + dset = h5_file.create_dataset('objects/'+object_name+'/spectrum', + data=data[:, 0:3], + dtype='f') + + dset.attrs['instrument'] = str(instrument) sys.stdout.write(' [DONE]\n') sys.stdout.flush() @@ -249,10 +300,14 @@ def add_object(self, def add_photometry(self, library): """ - :param library: Photometry library. - :type library: str + Parameters + ---------- + library : str + Photometry library. - :return: None + Returns + ------- + None """ h5_file = h5py.File(self.database, 'a') @@ -281,17 +336,21 @@ def add_calibration(self, """ Function for adding a calibration spectrum to the database. - :param filename: Filename with the calibration spectrum. The first column should contain - the wavelength (micron), the second column the flux density (W m-2 - micron-1), and the third column the error (W m-2 micron-1). - :type filename: str - :param tag: Tag name in the database. - :type tag: str - :param scaling: Scaling for the wavelength and flux as (scaling_wavelength, scaling_flux). - Not used if set to None. - :type scaling: tuple(float, float) + Parameters + ---------- + filename : str + Filename with the calibration spectrum. The first column should contain the wavelength + (micron), the second column the flux density (W m-2 micron-1), and the third column + the error (W m-2 micron-1). + tag : str + Tag name in the database. + scaling : tuple(float, float) + Scaling for the wavelength and flux as (scaling_wavelength, scaling_flux). Not used if + set to None. - :return: None + Returns + ------- + None """ if not scaling: @@ -330,10 +389,14 @@ def add_calibration(self, def add_spectrum(self, spectrum): """ - :param spectrum: Spectral library. - :type spectrum: str + Parameters + ---------- + spectrum : str + Spectral library. - :return: None + Returns + ------- + None """ h5_file = h5py.File(self.database, 'a') @@ -360,14 +423,18 @@ def add_votable(self, distance, filename): """ - :param object_name: Object name. - :type object_name: str - :param distance: Distance (pc). - :type distance: float - :param filename: Filename. - :type filename: str + Parameters + ---------- + object_name : str + Object name. + distance : float + Distance (pc). + filename : str + Filename. - :return: None + Returns + ------- + None """ # flux = {} @@ -443,22 +510,26 @@ def add_samples(self, modelpar, distance=None): """ - :param sampler: Ensemble sampler. - :type sampler: emcee.ensemble.EnsembleSampler - :param spectrum: Tuple with the spectrum type ('model' or 'calibration') and spectrum name - (e.g. 'drift-phoenix'). - :type spectrum: tuple(str, str) - :param tag: Database tag. - :type tag: str - :param chisquare: Maximum likelihood solution. Tuple with the chi-square value and related - parameter values. - :type chisquare: tuple(float, float) - :param modelpar: List with the model parameter names. - :type modelpar: list(str, ) - :param distance: Distance to the object (pc). Not used if set to None. - :type distance: float - - :return: None + Parameters + ---------- + sampler : emcee.ensemble.EnsembleSampler + Ensemble sampler. + spectrum : tuple(str, str) + Tuple with the spectrum type ('model' or 'calibration') and spectrum name (e.g. + 'drift-phoenix'). + tag : str + Database tag. + chisquare : tuple(float, float) + Maximum likelihood solution. Tuple with the chi-square value and related parameter + values. + modelpar : list(str, ) + List with the model parameter names. + distance : float + Distance to the object (pc). Not used if set to None. + + Returns + ------- + None """ h5_file = h5py.File(self.database, 'a') @@ -515,11 +586,12 @@ def add_samples(self, def get_chisquare(self, tag): """ - :param tag: - :type tag: str + Parameters + ---------- + tag : str - :return: - :rtype: species.core.box.SamplesBox + Returns + ------- """ h5_file = h5py.File(self.database, 'r') @@ -544,22 +616,30 @@ def get_mcmc_spectra(self, burnin, random, wavelength, - sampling=None): - """ - :param tag: Database tag with the MCMC samples. - :type tag: str - :param burnin: Number of burnin steps. - :type burnin: int - :param random: Number of random samples. - :type random: int - :param wavelength: Wavelength range (micron) or filter name. Full spectrum if set to None. - :type wavelength: tuple(float, float) or str - :param sampling: Spectral sampling. - :type sampling: tuple - - :return: Boxes with the randomly sampled spectra. - :rtype: tuple(species.core.box.ModelBox, ) - """ + specres=None): + """ + Parameters + ---------- + tag : str + Database tag with the MCMC samples. + burnin : int + Number of burnin steps. + random : int + Number of random samples. + wavelength : tuple(float, float) or str + Wavelength range (micron) or filter name. Full spectrum if set to None. + specres : float + Spectral resolution, achieved by smoothing with a Gaussian kernel. The original + wavelength points are used if set to None. + + Returns + ------- + tuple(species.core.box.ModelBox, ) + Boxes with the randomly sampled spectra. + """ + + sys.stdout.write('Getting MCMC spectra...') + sys.stdout.flush() h5_file = h5py.File(self.database, 'r') dset = h5_file['results/mcmc/'+tag] @@ -591,6 +671,10 @@ def get_mcmc_spectra(self, boxes = [] + progbar = progress.bar.Bar('\rGetting MCMC spectra...', + max=samples.shape[0], + suffix='%(percent)d%%') + for i in range(samples.shape[0]): model_par = {} for j in range(samples.shape[1]): @@ -600,7 +684,7 @@ def get_mcmc_spectra(self, model_par['distance'] = distance if spectrum_type == 'model': - specbox = readmodel.get_model(model_par, sampling) + specbox = readmodel.get_model(model_par, specres) elif spectrum_type == 'calibration': specbox = readcalib.get_spectrum(model_par) @@ -608,22 +692,30 @@ def get_mcmc_spectra(self, boxes.append(specbox) + progbar.next() + + progbar.finish() + h5_file.close() return tuple(boxes) def get_object(self, object_name, - filter_id): + filter_id=None): """ - :param object_name: Object name in the database. - :type object_name: str - :param filter_id: Filter IDs for which the photometry is selected. All available photometry - of the object is selected if set to None. - :type filter_id: tuple(str, ) + Parameters + ---------- + object_name : str + Object name in the database. + filter_id : tuple(str, ) + Filter IDs for which the photometry is selected. All available photometry of the object + is selected if set to None. - :return: Box with the object. - :rtype: species.core.box.ObjectBox + Returns + ------- + species.core.box.ObjectBox + Box with the object. """ h5_file = h5py.File(self.database, 'r') @@ -643,13 +735,18 @@ def get_object(self, else: for key in dset.keys(): - if key != 'distance': + 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]) + if 'objects/'+object_name+'/spectrum' in h5_file: + spectrum = np.asarray(h5_file['objects/'+object_name+'/spectrum']) + else: + spectrum = None + h5_file.close() return box.create_box('object', @@ -657,22 +754,24 @@ def get_object(self, filter=tuple(magnitude.keys()), magnitude=magnitude, flux=flux, - distance=distance) + distance=distance, + spectrum=spectrum) def get_samples(self, tag, burnin=None, random=None): """ - :param tag: - :type tag: str - :param burnin: - :type burnin: int - :param random: - :type random: int + Parameters + ---------- + tag: str + burnin : int + random : int - :return: - :rtype: species.core.box.SamplesBox + Returns + ------- + species.core.box.SamplesBox + Box with the MCMC samples. """ h5_file = h5py.File(self.database, 'r') diff --git a/species/plot/plot_mcmc.py b/species/plot/plot_mcmc.py index 4fca9674..a24c7137 100644 --- a/species/plot/plot_mcmc.py +++ b/species/plot/plot_mcmc.py @@ -1,5 +1,5 @@ """ -Module with functions for making plots. +Module for plotting MCMC results. """ import os @@ -25,7 +25,22 @@ def plot_walkers(tag, nsteps=None, offset=None): """ - :return: None + Function to plot the step history of the walkers. + + Parameters + ---------- + tag : str + Database tag with the MCMC samples. + output : str + Output filename. + nsteps : int + Number of steps. + offset : tuple(float, float) + Offset of the x- and y-axis label. + + Returns + ------- + None """ sys.stdout.write('Plotting walkers: '+output+'...') @@ -100,7 +115,26 @@ def plot_posterior(tag, offset=None, title_fmt='.2f'): """ - :return: None + Function to plot the posterior distributions. + + Parameters + ---------- + tag : str + Database tag with the MCMC samples. + burnin : int + Number of burnin steps to exclude. + output : str + Output filename. + title : str + Plot title. + offset : tuple(float, float) + Offset of the x- and y-axis label. + title_fmt : str + Format of the median and error values. + + Returns + ------- + None """ sys.stdout.write('Plotting posteriors: '+output+'...') diff --git a/species/plot/plot_spectrum.py b/species/plot/plot_spectrum.py index 90e620d9..ed8d2ce3 100644 --- a/species/plot/plot_spectrum.py +++ b/species/plot/plot_spectrum.py @@ -1,5 +1,5 @@ """ -Module with functions for making plots. +Module with a function for plotting spectra. """ import os @@ -32,24 +32,47 @@ def plot_spectrum(boxes, scale=('linear', 'linear'), title=None, offset=None, - legend='upper left'): + legend='upper left', + figsize=(7., 5.)): """ - :param boxes: - :type boxes: tuple(species.analysis.box.SpectrumBox and/or - species.analysis.box.PhotometryBox and/or - species.analysis.box.ModelBox) - :param filters: - :type filters: tuple(str, ) - :param output: - :type output: str - - :return: None + Parameters + ---------- + boxes : tuple(species.core.box, ) + Boxes with data. + filters : tuple(str, ) + Filter IDs for which the transmission profile is plotted. + output : str + Output filename. + colors : tuple(str, ) + 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. + residuals : species.core.box.ResidualsBox + Box with residuals of a fit. + xlim : tuple(float, float) + Limits of the x-axis. + ylim : tuple(float, float) + Limits of the y-axis. + scale : tuple(str, str) + Scale of the axes ('linear' or 'log'). + title : str + Title. + offset : tuple(float, float) + Offset for the label of the x- and y-axis. + legend : str + Location of the legend. + figsize : tuple(float, float) + Figure size. + + Returns + ------- + None """ marker = itertools.cycle(('o', 's', '*', 'p', '<', '>', 'P', 'v', '^')) if residuals and filters: - plt.figure(1, figsize=(7, 5)) + plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(3, 1, height_ratios=[1, 3, 1]) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) @@ -58,7 +81,7 @@ def plot_spectrum(boxes, ax3 = plt.subplot(gridsp[2, 0]) elif residuals: - plt.figure(1, figsize=(7, 5)) + plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(2, 1, height_ratios=[4, 1]) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) @@ -66,7 +89,7 @@ def plot_spectrum(boxes, ax3 = plt.subplot(gridsp[1, 0]) elif filters: - plt.figure(1, figsize=(7, 5)) + plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(2, 1, height_ratios=[1, 4]) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) @@ -74,7 +97,7 @@ def plot_spectrum(boxes, ax2 = plt.subplot(gridsp[0, 0]) else: - plt.figure(1, figsize=(7, 4)) + plt.figure(1, figsize=figsize) gridsp = mpl.gridspec.GridSpec(1, 1) gridsp.update(wspace=0, hspace=0, left=0, right=1, bottom=0, top=1) @@ -211,6 +234,9 @@ def plot_spectrum(boxes, if residuals: ax3.set_xscale(scale[0]) + color_obj_phot = None + color_obj_spec = None + for j, boxitem in enumerate(boxes): if isinstance(boxitem, (box.SpectrumBox, box.ModelBox)): wavelength = boxitem.wavelength @@ -247,7 +273,8 @@ def plot_spectrum(boxes, label = None if colors: - ax1.plot(wavelength, masked/scaling, color=colors[j], lw=0.5, label=label, zorder=4) + ax1.plot(wavelength, masked/scaling, color=colors[j], lw=0.5, + label=label, zorder=4) else: ax1.plot(wavelength, masked/scaling, lw=0.5, label=label, zorder=4) @@ -267,7 +294,8 @@ def plot_spectrum(boxes, 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=3) + ax1.plot(wavelength, masked/scaling, lw=0.2, color=colors[j], + alpha=0.5, zorder=3) else: ax1.plot(wavelength, masked/scaling, lw=0.2, alpha=0.5, zorder=3) @@ -285,9 +313,25 @@ def plot_spectrum(boxes, wavelength = transmission.mean_wavelength() fwhm = transmission.filter_fwhm() + color_obj_phot = colors[j][0] + ax1.errorbar(wavelength, boxitem.flux[item][0]/scaling, xerr=fwhm/2., - yerr=boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=5, - color=colors[j], markerfacecolor=colors[j]) + yerr=boxitem.flux[item][1]/scaling, marker='s', ms=5, zorder=6, + color=color_obj_phot, markerfacecolor=color_obj_phot) + + if boxitem.spectrum is not None: + masked = np.ma.array(boxitem.spectrum, mask=np.isnan(boxitem.spectrum)) + + color_obj_spec = colors[j][1] + + if colors is None: + ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling, + ms=2, marker='s', zorder=5) + + else: + ax1.errorbar(masked[:, 0], masked[:, 1]/scaling, yerr=masked[:, 2]/scaling, + marker='o', ms=2, zorder=5, color=color_obj_spec, + markerfacecolor=color_obj_spec) elif isinstance(boxitem, box.SynphotBox): for item in boxitem.flux: @@ -296,7 +340,7 @@ def plot_spectrum(boxes, fwhm = transmission.filter_fwhm() ax1.errorbar(wavelength, boxitem.flux[item]/scaling, xerr=fwhm/2., yerr=None, - alpha=0.7, marker='s', ms=5, zorder=6, color=colors[j], + alpha=0.7, marker='s', ms=5, zorder=7, color=color_obj_spec, markerfacecolor='white') if filters: @@ -307,17 +351,24 @@ def plot_spectrum(boxes, ax2.plot(data[0, ], data[1, ], '-', lw=0.7, color='black') if residuals: - diff = np.zeros((2, len(residuals[0].flux))) + res_max = 0. - for i, item in enumerate(residuals[0].flux): - transmission = read_filter.ReadFilter(item) - diff[0, i] = transmission.mean_wavelength() - diff[1, i] = (residuals[0].flux[item][0]-residuals[1].flux[item]) / \ - residuals[0].flux[item][1] + 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) + + res_max = np.nanmax(np.abs(residuals.photometry[1, ])) - ax3.plot(diff[0, ], diff[1, ], marker='s', ms=5, linestyle='none', color='black') + if residuals.spectrum is not None: + ax3.plot(residuals.spectrum[0, ], residuals.spectrum[1, ], marker='o', + ms=2, linestyle='none', color=color_obj_spec, zorder=1) - res_lim = math.ceil(np.amax(np.abs(diff[1, ]))) + max_tmp = np.nanmax(np.abs(residuals.spectrum[1, ])) + + if max_tmp > res_max: + res_max = max_tmp + + res_lim = math.ceil(max_tmp) ax3.axhline(0.0, linestyle='--', color='gray', dashes=(2, 4), zorder=1) ax3.set_ylim(-res_lim, res_lim) @@ -328,9 +379,6 @@ def plot_spectrum(boxes, sys.stdout.write('Plotting spectrum: '+output+'...') sys.stdout.flush() - # ax1.xaxis.set_major_formatter(FormatStrFormatter('%.1f')) - # ax1.xaxis.set_minor_formatter(FormatStrFormatter('%.1f')) - if title: if filters: ax2.set_title(title, y=1.02, fontsize=15) diff --git a/species/read/read_model.py b/species/read/read_model.py index 037c7251..3cdfab0a 100644 --- a/species/read/read_model.py +++ b/species/read/read_model.py @@ -1,5 +1,5 @@ """ -Text +Module for reading atmospheric models. """ import os @@ -19,7 +19,7 @@ class ReadModel: """ - Text + Read atmospheric model spectra. """ def __init__(self, @@ -32,7 +32,7 @@ def __init__(self, model : str Model name. wavelength : tuple(float, float) or str - Wavelength range (micron) or filter name. Full spectrum if set to None. + Wavelength range (micron) or filter name. Full spectrum is used if set to None. teff : tuple(float, float) Effective temperature (K) range. Restricting the temperature range will speed up the computation. @@ -42,10 +42,12 @@ def __init__(self, None """ - self.model = model self.teff = teff + self.spectrum_interp = None + self.wl_points = None + self.wl_index = None if isinstance(wavelength, str): self.filter_name = wavelength @@ -56,13 +58,6 @@ def __init__(self, self.filter_name = None self.wavelength = wavelength - config_file = os.path.join(os.getcwd(), 'species_config.ini') - - config = configparser.ConfigParser() - config.read_file(open(config_file)) - - self.database = config['species']['database'] - def open_database(self): """ Returns @@ -71,7 +66,14 @@ def open_database(self): Database. """ - h5_file = h5py.File(self.database, 'r') + config_file = os.path.join(os.getcwd(), 'species_config.ini') + + config = configparser.ConfigParser() + config.read_file(open(config_file)) + + database_path = config['species']['database'] + + h5_file = h5py.File(database_path, 'r') try: h5_file['models/'+self.model] @@ -80,10 +82,44 @@ def open_database(self): h5_file.close() species_db = database.Database() species_db.add_model(self.model, self.wavelength, self.teff) - h5_file = h5py.File(self.database, 'r') + h5_file = h5py.File(database_path, 'r') return h5_file + def wavelength_points(self, + hdf5_file): + """ + Parameters + ---------- + hdf5_file : h5py._hl.files.File + hdf5_file. + + Returns + ------- + numpy.ndarray + Wavelength points (micron). + numpy.ndarray + Array with the size of the original wavelength grid. The booleans indicate if a + wavelength point was used. + """ + + wl_points = np.asarray(hdf5_file['models/'+self.model+'/wavelength']) + + if self.wavelength is None: + wl_index = np.ones(wl_points.shape[0], dtype=bool) + + else: + wl_index = (wl_points > self.wavelength[0]) & (wl_points < self.wavelength[1]) + index = np.where(wl_index)[0] + + if index[0]-1 >= 0: + wl_index[index[0] - 1] = True + + if index[-1]+1 < wl_index.size: + wl_index[index[-1] + 1] = True + + return wl_points[wl_index], wl_index + def interpolate(self): """ Returns @@ -93,49 +129,40 @@ def interpolate(self): h5_file = self.open_database() - wavelength = np.asarray(h5_file['models/'+self.model+'/wavelength']) flux = np.asarray(h5_file['models/'+self.model+'/flux']) teff = np.asarray(h5_file['models/'+self.model+'/teff']) logg = np.asarray(h5_file['models/'+self.model+'/logg']) - wl_index = (wavelength > self.wavelength[0]) & (wavelength < self.wavelength[1]) + if self.wl_points is None: + self.wl_points, self.wl_index = self.wavelength_points(h5_file) + + if self.model in ('drift-phoenix', 'bt-nextgen', 'petitcode_warm_clear'): + feh = np.asarray(h5_file['models/'+self.model+'/feh']) - index = np.where(wl_index)[0] + points = (teff, logg, feh) + flux = flux[:, :, :, self.wl_index] - if index[0]-1 >= 0: - wl_index[index[0] - 1] = True + elif self.model == 'petitcode_warm_cloudy': + feh = np.asarray(h5_file['models/petitcode_warm_cloudy/feh']) + fsed = np.asarray(h5_file['models/petitcode_warm_cloudy/fsed']) - if index[-1]+1 < wl_index.size: - wl_index[index[-1] + 1] = True + points = (teff, logg, feh, fsed) + flux = flux[:, :, :, :, self.wl_index] - if self.model in ('drift-phoenix', 'bt-nextgen'): - feh = np.asarray(h5_file['models/'+self.model+'/feh']) - flux = flux[:, :, :, wl_index] - points = (teff, logg, feh, wavelength[wl_index]) - - # elif self.model == 'petitcode_warm_clear': - # feh = np.asarray(h5_file['models/petitcode_warm_clear/feh']) - # flux = flux[:, :, :, wl_index] - # points = np.asarray((teff, logg, feh, wavelength[wl_index])) - # - # elif self.model == 'petitcode_warm_cloudy': - # feh = np.asarray(h5_file['models/petitcode_warm_cloudy/feh']) - # fsed = np.asarray(h5_file['models/petitcode_warm_cloudy/fsed']) - # flux = flux[:, :, :, :, wl_index] - # points = np.asarray((teff, logg, feh, fsed, wavelength[wl_index])) - # - # elif self.model == 'petitcode_hot_clear': - # feh = np.asarray(h5_file['models/petitcode_hot_clear/feh']) - # co_ratio = np.asarray(h5_file['models/petitcode_hot_clear/co']) - # flux = flux[:, :, :, :, wl_index] - # points = np.asarray((teff, logg, feh, co_ratio, wavelength[wl_index])) - # - # elif self.model == 'petitcode_hot_cloudy': - # feh = np.asarray(h5_file['models/petitcode_hot_cloudy/feh']) - # co_ratio = np.asarray(h5_file['models/petitcode_hot_cloudy/co']) - # fsed = np.asarray(h5_file['models/petitcode_hot_cloudy/fsed']) - # flux = flux[:, :, :, :, :, wl_index] - # points = np.asarray((teff, logg, feh, co_ratio, fsed, wavelength[wl_index])) + elif self.model == 'petitcode_hot_clear': + feh = np.asarray(h5_file['models/petitcode_hot_clear/feh']) + co_ratio = np.asarray(h5_file['models/petitcode_hot_clear/co']) + + points = (teff, logg, feh, co_ratio) + flux = flux[:, :, :, :, self.wl_index] + + elif self.model == 'petitcode_hot_cloudy': + feh = np.asarray(h5_file['models/petitcode_hot_cloudy/feh']) + co_ratio = np.asarray(h5_file['models/petitcode_hot_cloudy/co']) + fsed = np.asarray(h5_file['models/petitcode_hot_cloudy/fsed']) + + flux = flux[:, :, :, :, :, self.wl_index] + points = (teff, logg, feh, co_ratio, fsed) h5_file.close() @@ -145,6 +172,83 @@ def interpolate(self): bounds_error=False, fill_value=np.nan) + def get_model(self, + model_par, + specres=None): + """ + Parameters + ---------- + model_par : dict + Model parameter values. + specres : float + Spectral resolution, achieved by smoothing with a Gaussian kernel. The original + wavelength points are used if set to None. + + Returns + ------- + species.core.box.ModelBox + Box with the model spectrum. + """ + + if self.spectrum_interp is None: + self.interpolate() + + if self.wavelength is None: + wl_points = self.get_wavelength() + self.wavelength = (wl_points[0], wl_points[-1]) + + if self.model in ('drift-phoenix', 'bt-nextgen', 'petitcode_warm_clear'): + parameters = [model_par['teff'], + model_par['logg'], + model_par['feh']] + + elif self.model == 'petitcode_warm_cloudy': + parameters = [model_par['teff'], + model_par['logg'], + model_par['feh'], + model_par['fsed']] + + elif self.model == 'petitcode_hot_clear': + parameters = [model_par['teff'], + model_par['logg'], + model_par['feh'], + model_par['co']] + + elif self.model == 'petitcode_hot_cloudy': + parameters = [model_par['teff'], + model_par['logg'], + model_par['feh'], + model_par['co'], + model_par['fsed']] + + flux = self.spectrum_interp(parameters)[0] + + if 'radius' in model_par: + model_par['mass'] = read_util.get_mass(model_par) + + if 'distance' in model_par: + scaling = (model_par['radius']*constants.R_JUP)**2 / \ + (model_par['distance']*constants.PARSEC)**2 + + flux *= scaling + + if specres is not None: + index = np.where(np.isnan(flux))[0] + + if index.size > 0: + raise ValueError('Flux values should not contains NaNs.') + + flux = read_util.smooth_spectrum(wavelength=self.wl_points, + flux=flux, + specres=specres, + size=11) + + return box.create_box(boxtype='model', + model=self.model, + wavelength=self.wl_points, + flux=flux, + parameters=model_par) + def get_data(self, model_par): """ @@ -157,26 +261,17 @@ def get_data(self, Returns ------- species.core.box.ModelBox - Spectrum (micron, W m-2 micron-1). + Box with the model spectrum. """ h5_file = self.open_database() - wavelength = np.asarray(h5_file['models/'+self.model+'/wavelength']) + wl_points, wl_index = self.wavelength_points(h5_file) + flux = np.asarray(h5_file['models/'+self.model+'/flux']) teff = np.asarray(h5_file['models/'+self.model+'/teff']) logg = np.asarray(h5_file['models/'+self.model+'/logg']) - if self.wavelength is None: - wl_index = np.ones(wavelength.shape[0], dtype=bool) - - else: - wl_index = (wavelength > self.wavelength[0]) & (wavelength < self.wavelength[1]) - index = np.where(wl_index)[0] - - wl_index[index[0] - 1] = True - wl_index[index[-1] + 1] = True - if self.model in ('drift-phoenix', 'bt-nextgen'): feh = np.asarray(h5_file['models/'+self.model+'/feh']) @@ -201,182 +296,83 @@ def get_data(self, else: feh_index = feh_index[0] - wavelength = wavelength[wl_index] flux = flux[teff_index, logg_index, feh_index, wl_index] if 'radius' in model_par and 'distance' in model_par: scaling = (model_par['radius']*constants.R_JUP)**2 / \ (model_par['distance']*constants.PARSEC)**2 + flux *= scaling + h5_file.close() + return box.create_box(boxtype='model', model=self.model, - wavelength=wavelength, + wavelength=wl_points, flux=flux, parameters=model_par) - def get_model(self, - model_par, - sampling): + def get_photometry(self, + model_par, + synphot=None): """ Parameters ---------- model_par : dict Model parameter values. - sampling : tuple - Type of wavelength sampling. + synphot : species.analysis.photometry.SyntheticPhotometry + Synthetic photometry object. Returns ------- - species.core.box.ModelBox - Spectrum (micron, W m-2 micron-1). + float + Average flux density (W m-2 micron-1). """ - if not self.wavelength: - wl_points = self.get_wavelength() - self.wavelength = (wl_points[0], wl_points[-1]) - if self.spectrum_interp is None: self.interpolate() - wavelength = [self.wavelength[0]] - - if sampling[0] == 'specres': - while wavelength[-1] <= self.wavelength[1]: - wavelength.append(wavelength[-1] + wavelength[-1]/sampling[1]) - wavelength = np.asarray(wavelength[:-1]) - - elif sampling[0] == 'gaussian': - wavelength = np.linspace(self.wavelength[0], self.wavelength[1], sampling[1][0]) - - flux = np.zeros(wavelength.shape) - - for i, item in enumerate(wavelength): - - if self.model in ('drift-phoenix', 'bt-nextgen'): - parameters = [model_par['teff'], - model_par['logg'], - model_par['feh'], - item] - - # elif self.model == 'petitcode_warm_clear': - # parameters = [model_par['teff'], - # model_par['logg'], - # model_par['feh'], - # item] - # - # elif self.model == 'petitcode_warm_cloudy': - # parameters = [model_par['teff'], - # model_par['logg'], - # model_par['feh'], - # model_par['fsed'], - # item] - # - # elif self.model == 'petitcode_hot_clear': - # parameters = [model_par['teff'], - # model_par['logg'], - # model_par['feh'], - # model_par['co'], - # item] - # - # elif self.model == 'petitcode_hot_cloudy': - # parameters = [model_par['teff'], - # model_par['logg'], - # model_par['feh'], - # model_par['co'], - # model_par['fsed'], - # item] - - flux[i] = self.spectrum_interp(np.asarray(parameters)) - - if 'radius' in model_par: - model_par['mass'] = read_util.get_mass(model_par) - - if 'distance' in model_par: - scaling = (model_par['radius']*constants.R_JUP)**2 / \ - (model_par['distance']*constants.PARSEC)**2 - - flux *= scaling - - if sampling[0] == 'gaussian': - index = np.where(np.isnan(flux))[0] + spectrum = self.get_model(model_par, None) - wavelength = np.delete(wavelength, index) - flux = np.delete(flux, index) - - flux = read_util.smooth_spectrum(wavelength, flux, sampling[1][1], 11) + if not synphot: + synphot = photometry.SyntheticPhotometry(self.filter_name) - return box.create_box(boxtype='model', - model=self.model, - wavelength=wavelength, - flux=flux, - parameters=model_par) + return synphot.spectrum_to_photometry(spectrum.wavelength, spectrum.flux) - def get_photometry(self, - model_par, - sampling, - synphot=None): + def get_magnitude(self, + model_par): """ Parameters ---------- model_par : dict Model parameter values. - sampling : float - Spectral sampling. The original grid is used (nearest model parameter values) if set - to none. - synphot : species.analysis.photometry.SyntheticPhotometry - Synthetic photometry object. Returns ------- float - Average flux density (W m-2 micron-1). + Apparent magnitude (mag). + float + Absolute magnitude (mag). """ - if sampling is None: - spectrum = self.get_data(model_par) + if self.spectrum_interp is None: + self.interpolate() - else: - if self.spectrum_interp is None: - self.interpolate() + spectrum = self.get_model(model_par, None) - spectrum = self.get_model(model_par, sampling) + synphot = photometry.SyntheticPhotometry(self.filter_name) - if not synphot: - synphot = photometry.SyntheticPhotometry(self.filter_name) + if 'distance' in model_par: + app_mag, abs_mag = synphot.spectrum_to_magnitude(spectrum.wavelength, + spectrum.flux, + model_par['distance']) - return synphot.spectrum_to_photometry(spectrum.wavelength, spectrum.flux) + else: + app_mag, abs_mag = synphot.spectrum_to_magnitude(spectrum.wavelength, + spectrum.flux, + None) - # def get_magnitude(self, model_par, sampling): - # """ - # :param model_par: Model parameter values. - # :type model_par: dict - # :param sampling: Spectral sampling. The original grid is used (nearest model parameter - # values) if set to none. - # :type sampling: float - # - # :return: Apparent magnitude (mag), absolute magnitude (mag). - # :rtype: float, float - # """ - # - # if sampling is None: - # spectrum = self.get_data(model_par) - # - # else: - # if self.spectrum_interp is None: - # self.interpolate() - # - # spectrum = self.get_model(model_par, sampling) - # - # transmission = read_filter.ReadFilter(self.filter_name) - # filter_interp = transmission.interpolate() - # - # synphot = photometry.SyntheticPhotometry(filter_interp) - # mag = synphot.spectrum_to_magnitude(spectrum.wavelength, - # spectrum.flux, - # model_par['distance']) - # - # return mag[0], mag[1] + return app_mag, abs_mag def get_bounds(self): """ @@ -458,4 +454,6 @@ def get_parameters(self): for i in range(nparam): param.append(dset.attrs['parameter'+str(i)]) + h5_file.close() + return param diff --git a/species/read/read_object.py b/species/read/read_object.py index afefa8ee..e2f1f808 100644 --- a/species/read/read_object.py +++ b/species/read/read_object.py @@ -1,5 +1,5 @@ """ -Read module. +Module for reading object data. """ import os @@ -13,16 +13,20 @@ class ReadObject: """ - Text + Reading object data from the database. """ def __init__(self, object_name): """ - :param object_name: Object name. - :type object_name: str - - :return: None + Parameters + ---------- + object_name : str + Object name. + + Returns + ------- + None """ self.object_name = object_name @@ -37,12 +41,16 @@ def __init__(self, def get_photometry(self, filter_name): """ - :param filter_name: Filter name. - :type filter_name: str - - :return: Apparent magnitude (mag), magnitude error (error), - apparent flux (W m-2 micron-1), flux error (W m-2 micron-1). - :rtype: numpy.ndarray + Parameters + ---------- + filter_name : str + Filter name. + + Returns + ------- + numpy.ndarray + Apparent magnitude (mag), magnitude error (error), apparent flux (W m-2 micron-1), + flux error (W m-2 micron-1). """ h5_file = h5py.File(self.database, 'r') @@ -51,10 +59,41 @@ def get_photometry(self, return obj_phot + def get_spectrum(self): + """ + Returns + ------- + numpy.ndarray + Wavelength (micron), apparent flux (W m-2 micron-1), and flux error (W m-2 micron-1). + """ + + h5_file = h5py.File(self.database, 'r') + spectrum = np.asarray(h5_file['objects/'+self.object_name+'/spectrum']) + h5_file.close() + + return spectrum + + def get_instrument(self): + """ + Returns + ------- + str + Instrument that was used for the spectrum. + """ + + h5_file = h5py.File(self.database, 'r') + dset = h5_file['objects/'+self.object_name+'/spectrum'] + instrument = dset.attrs['instrument'] + h5_file.close() + + return instrument + def get_distance(self): """ - :return: Distance (pc). - :rtype: float + Returns + ------- + float + Distance (pc). """ h5_file = h5py.File(self.database, 'r') @@ -66,11 +105,15 @@ def get_distance(self): def get_absmag(self, filter_name): """ - :param filter_name: Filter name. - :type filter_name: str - - :return: Absolute magnitude (mag), magnitude error (error). - :rtype: float, float + Parameters + ---------- + filter_name : str + Filter name. + + Returns + ------- + float, float + Absolute magnitude (mag), magnitude error (error). """ h5_file = h5py.File(self.database, 'r') diff --git a/species/util/data_util.py b/species/util/data_util.py index e10bed2e..5d24f3e5 100644 --- a/species/util/data_util.py +++ b/species/util/data_util.py @@ -1,5 +1,5 @@ """ -Utility functions. +Utility functions for data processing. """ import numpy as np @@ -9,11 +9,15 @@ def update_sptype(sptypes): """ Function to update a list with spectral types to two characters (e.g., M8, L3, or T1). - Args: - sptypes(numpy.ndarray): Spectral types. + Parameters + ---------- + sptypes : numpy.ndarray + Spectral types. - Returns: - numpy.ndarray: Updated spectral types. + Returns + ------- + numpy.ndarray + Updated spectral types. """ mlty = ('O', 'B', 'A', 'F', 'G', 'K', 'M', 'L', 'T', 'Y') @@ -42,11 +46,15 @@ def update_filter(filter_in): Function to update filter ID from the Vizier Photometry viewer VOTable to the filter ID from the SVO Filter Profile Service. - Args: - filter_in(str): Filter ID in the format of the Vizier Photometry viewer + Parameters + ---------- + filter_in : str + Filter ID in the format of the Vizier Photometry viewer. - Returns: - str: Filter ID in the format of the SVO Filter Profile Service. + Returns + ------- + str + Filter ID in the format of the SVO Filter Profile Service. """ if filter_in[0:5] == b'2MASS': @@ -70,15 +78,21 @@ def sort_data(teff, wavelength, flux): """ - Args: - teff(numpy.ndarray) - logg(numpy.ndarray) - feh(numpy.ndarray) - wavelength(numpy.ndarray) - flux(numpy.ndarray) - - Returns: - tuple(numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray) + Parameters + ---------- + teff : numpy.ndarray + logg : numpy.ndarray + feh : numpy.ndarray + wavelength : numpy.ndarray + flux : numpy.ndarray + + Returns + ------- + numpy.ndarray + numpy.ndarray + numpy.ndarray + numpy.ndarray + numpy.ndarray """ teff_unique = np.unique(teff) @@ -104,12 +118,17 @@ def write_data(model, database, data_sorted): """ - Args: - model(str): Atmosphere model. - database(h5py._hl.files.File) - - Returns: - None + Parameters + ---------- + model : str + Atmosphere model. + database: h5py._hl.files.File + Database. + data_sorted : + + Returns + ------- + None """ if 'models/'+model in database: @@ -146,12 +165,16 @@ def write_data(model, def add_missing(model, database): """ - Args: - model(str): Atmosphere model. - database(h5py._hl.files.File) - - Returns: - None + Parameters + ---------- + model : str + Atmosphere model. + database : h5py._hl.files.File + Database. + + Returns + ------- + None """ teff = np.asarray(database['models/'+model+'/teff']) diff --git a/species/util/phot_util.py b/species/util/phot_util.py new file mode 100644 index 00000000..05fa9d0c --- /dev/null +++ b/species/util/phot_util.py @@ -0,0 +1,133 @@ +""" +Utility functions for photometry. +""" + +import spectres +import numpy as np + +from species.core import box +from species.read import read_model, read_calibration, read_filter + + +def multi_photometry(datatype, + spectrum, + filters, + model_par): + """ + Parameters + ---------- + datatype : str + Data type ('model' or 'calibration'). + spectrum : str + Spectrum name (e.g., 'drift-phoenix'). + filters : tuple(str, ) + Filter IDs. + model_par : dict + Model parameter values. + + Returns + ------- + species.core.box.SynphotBox + Box with synthetic photometry. + """ + + flux = {} + + if datatype == 'model': + for item in filters: + readmodel = read_model.ReadModel(spectrum, item) + flux[item] = readmodel.get_photometry(model_par) + + elif datatype == 'calibration': + for item in filters: + readcalib = read_calibration.ReadCalibration(spectrum, item) + flux[item] = readcalib.get_photometry(model_par) + + return box.create_box('synphot', name='synphot', flux=flux) + + +def apparent_to_absolute(app_mag, + distance): + """ + Parameters + ---------- + app_mag : float or numpy.ndarray + Apparent magnitude (mag). + distance : float or numpy.ndarray + Distance (pc). + + Returns + ------- + float or numpy.ndarray + Absolute magnitude (mag). + """ + + return app_mag - 5.*np.log10(distance) + 5. + +def get_residuals(model, + model_par, + filters, + objectbox): + """ + Parameters + ---------- + model : str + Atmospheric model. + model_par : dict + Model parameters and values. + filters : tuple(str, ) + Filter IDs. All available photometry of the object is used if set to None. + objectbox : species.core.box.ObjectBox + Box with the photometry and/or spectrum of an object. + + Returns + ------- + species.core.box.ResidualsBox + Box with the photometry and/or spectrum residuals. + """ + + if filters is None: + filters = objectbox.filter + + model_phot = multi_photometry(datatype='model', + spectrum=model, + filters=filters, + model_par=model_par) + + if objectbox.flux is not None: + res_phot = np.zeros((2, len(objectbox.flux))) + + for i, item in enumerate(objectbox.flux): + transmission = read_filter.ReadFilter(item) + + res_phot[0, i] = transmission.mean_wavelength() + res_phot[1, i] = (objectbox.flux[item][0]-model_phot.flux[item])/objectbox.flux[item][1] + + else: + res_phot = None + + if objectbox.spectrum is not None: + wl_range = (0.9*objectbox.spectrum[0, 0], 1.1*objectbox.spectrum[-1, 0]) + + readmodel = read_model.ReadModel(model, wl_range) + model = readmodel.get_model(model_par) + + wl_new = objectbox.spectrum[:, 0] + + flux_new = spectres.spectres(new_spec_wavs=wl_new, + old_spec_wavs=model.wavelength, + spec_fluxes=model.flux, + spec_errs=None) + + res_spec = np.zeros((2, objectbox.spectrum.shape[0])) + + res_spec[0, :] = wl_new + res_spec[1, :] = (objectbox.spectrum[:, 1]-flux_new)/objectbox.spectrum[:, 2] + + else: + res_spec = None + + return box.create_box(boxtype='residuals', + name=objectbox.name, + photometry=res_phot, + spectrum=res_spec) diff --git a/species/util/plot_util.py b/species/util/plot_util.py index a35a2a75..3c8e8f67 100644 --- a/species/util/plot_util.py +++ b/species/util/plot_util.py @@ -1,5 +1,5 @@ """ -Text. +Utility functions for plotting data. """ import numpy as np diff --git a/species/util/read_util.py b/species/util/read_util.py index aec391a4..70fc6f1c 100644 --- a/species/util/read_util.py +++ b/species/util/read_util.py @@ -57,7 +57,7 @@ def add_luminosity(modelbox, """ readmodel = read_model.ReadModel(model=modelbox.model, wavelength=None, teff=None) - fullspec = readmodel.get_model(model_par=modelbox.parameters, sampling=('specres', specres)) + fullspec = readmodel.get_model(model_par=modelbox.parameters) flux = simps(fullspec.flux, fullspec.wavelength)