Source code for sbpy.calib.core

# Licensed under a 3-clause BSD style license - see LICENSE.rst

__all__ = [
    'solar_spectrum',
    'vega_spectrum',
    'solar_fluxd',
    'vega_fluxd',
    'SpectralStandard',
    'Sun',
    'Vega',
    'FilterLookupError',
    'UndefinedSourceError'
]

__doctest_requires__ = {
    'SpectralStandard': ['synphot'],
    'solar_spectrum': ['synphot'],
    'vega_spectrum': ['synphot'],
    'Sun': ['synphot'],
    'Vega': ['synphot'],
}

import os
from abc import ABC
from warnings import warn
from functools import wraps
import inspect
import json

import numpy as np
from astropy.utils.state import ScienceState
from astropy.utils.data import get_pkg_data_filename
from astropy.table import Table, QTable
from astropy.io import ascii
import astropy.units as u

from ..spectroscopy.sources import SpectralSource, SynphotRequired
from ..exceptions import SbpyException, OptionalPackageUnavailable
from .. import bib
from . import solar_sources, vega_sources

try:
    import synphot
    from synphot import SpectralElement
except ImportError:
    synphot = None

    class SpectralElement:
        pass


[docs]class UndefinedSourceError(SbpyException): "SpectralStandard was initialized without a source, but it was accessed."
[docs]class FilterLookupError(SbpyException): "Attempted to look up filter in, e.g., solar_fluxd, but not present."
[docs]class SpectralStandard(SpectralSource, ABC): """Abstract base class for SBPy spectral standards. Parameters ---------- spectrum_state : `~astropy.utils.state.ScienceState` Context manager for this source's spectrum. fluxd_state : `~astropy.utils.state.ScienceState` Context manager for this source's calibration by filter. source : `~synphot.SourceSpectrum` or ``None`` The source spectrum or ``None`` if unspecified. description : string, optional A brief description of the source spectrum. bibcode : string, optional Bibliography code for citation (see `sbpy.bib.register`). Attributes ---------- wave - Wavelengths of the source spectrum. fluxd - Source spectrum. description - Brief description of the source spectrum. meta - Meta data from `source`, if any. """ def __init__(self, source, description=None, bibcode=None): self._source = source self._description = description self._bibcode = bibcode self._bibtask = '.'.join((self.__module__, self.__class__.__name__)) def __repr__(self): if self.description is None: return '<{}>'.format(self.__class__.__name__) else: return '<{}: {}>'.format(self.__class__.__name__, self.description)
[docs] @classmethod def from_builtin(cls, name): """Spectrum from a built-in `sbpy` source. Parameters ---------- name : string The name of the spectrum. See `show_builtin()` for available sources. """ from astropy.utils.data import _is_url try: parameters = getattr(cls._sources, name).copy() except AttributeError: msg = 'Unknown spectrum "{}". Valid spectra:\n'.format( name) + cls.show_builtin(print=False) raise UndefinedSourceError(msg) if not _is_url(parameters['filename']): # find in the module's location parameters['filename'] = get_pkg_data_filename( os.path.join('data', parameters['filename'])) return cls.from_file(**parameters)
[docs] @classmethod def from_default(cls): """Initialize new spectral standard from current default. The spectrum will be ``None`` if `synphot` is not available. """ if synphot is None: warn(OptionalPackageUnavailable( 'synphot is not installed, returning an empty spectral' ' standard.')) standard = cls(None) else: standard = cls._spectrum_state.get() return standard
[docs] @classmethod def show_builtin(cls, print=True): """List built-in spectra.""" from builtins import print as print_func sources = inspect.getmembers( cls._sources, lambda v: isinstance(v, dict)) rows = [] for name, source in sources: rows.append((name, source['description'])) tab = Table(rows=rows, names=('name', 'description')).pformat( max_lines=-1, max_width=-1) tab = '\n'.join(tab) if print: print_func(tab) else: return tab
@property def source(self): if self._source is None: raise UndefinedSourceError('The source is not defined.') if self._bibcode is not None: bib.register(self._bibtask, {self._description: self._bibcode}) return self._source
[docs] def observe(self, wfb, unit=None, interpolate=False, **kwargs): """Observe as through filters or spectrometer. Calls ``observe_bandpass``, ``observe_spectrum``, or ``self()``, as appropriate. Parameters ---------- wfb : `~astropy.units.Quantity`, `~synphot.SpectralElement`, string Wavelengths, frequencies, or bandpasses. Bandpasses may be a filter name (string). May also be a list of ``SpectralElement`` or strings. unit : string, `~astropy.units.Unit`, optional Units of the output (spectral flux density). interpolate : bool, optional For wavelengths/frequencies, set to ``True`` for interpolation instead of rebinning. Use this when the spectral resolution of the source is close to that of the requested wavelengths. **kwargs Additional keyword arguments for `~synphot.observation.Observation`, e.g., ``force``. Returns ------- fluxd : `~astropy.units.Quantity` """ if isinstance(wfb, (list, tuple)): fluxd = [] for i in range(len(wfb)): fluxd.append(self.observe(wfb[i], unit=unit, interpolate=interpolate, **kwargs)) fluxd = u.Quantity(fluxd) elif isinstance(wfb, str): lambda_eff, lambda_pivot, fluxd = self.observe_filter_name( wfb, unit=unit) elif isinstance(wfb, u.Quantity): if interpolate: fluxd = self(wfb, unit=unit) else: fluxd = self.observe_spectrum(wfb, unit=unit, **kwargs) elif isinstance(wfb, SpectralElement): lambda_eff, fluxd = self.observe_bandpass(wfb, unit=unit, **kwargs) else: raise TypeError('Unsupported type for `wfb` type: {}' .format(type(wfb))) return fluxd
@wraps(SpectralSource.observe_bandpass) def observe_bandpass(self, *args, **kwargs): if synphot is None: raise SynphotRequired( 'synphot is required for observations through bandpass') return super().observe_bandpass(*args, **kwargs) @wraps(SpectralSource.observe_spectrum) def observe_spectrum(self, *args, **kwargs): if synphot is None: raise SynphotRequired( 'synphot is required for observations through bandpass') return super().observe_spectrum(*args, **kwargs)
[docs] def observe_filter_name(self, filt, unit=None): """Flux density through this filter. Does not use the spectrum, but instead the flux density calibration manager. If the name of the filter is BP, then the expected keys are: BP : flux density BP(lambda eff) : effective wavelength, optional BP(lambda pivot) : pivot wavelength, optional Parameters ---------- filt : string Name of the filter. unit : string, `~astropy.units.Unit`, optional Spectral flux density units for the output. Returns ------- lambda_eff: `~astropy.units.Quantity` Effective wavelength. ``None`` if it is not provided. lambda_pivot: `~astropy.units.Quantity` Pivot wavelength. ``None`` if it is not provided. fluxd : `~astropy.units.Quantity` Spectral flux density. Raises ------ ``FilterLookupError`` if the filter is not defined. """ from .. import units as sbu source_fluxd = self._fluxd_state.get() try: fluxd = source_fluxd[filt] except KeyError: raise FilterLookupError( 'Filter "{}" is not present in `{}`' .format(filt, self._fluxd_state.__name__)) lambda_eff = source_fluxd.get(filt + '(lambda eff)') lambda_pivot = source_fluxd.get(filt + '(lambda pivot)') if unit is not None: unit = u.Unit(unit) # convert to requested units, may need lambda_pivot if lambda_pivot is None: equiv = [] else: equiv = u.spectral_density(lambda_pivot) # are VEGAmag involved and is a conversion needed? if (sbu.VEGA.is_equivalent((unit, fluxd.unit)) and not unit.is_equivalent(fluxd.unit)): equiv += sbu.spectral_density_vega(filt) try: fluxd = fluxd.to(unit, equiv) except u.UnitConversionError as e: raise type(e)( '{} Is "{}(lambda pivot)" required and' ' was it provided?'.format(e, filt)) return lambda_eff, lambda_pivot, fluxd
[docs] def color_index(self, wfb, unit): """Color index (magnitudes) and effective wavelengths. Parameters ---------- wfb : `~astropy.units.Quantity`, or tuple/list of `~synphot.SectralElement`, string Two wavelengths, frequencies, or bandpasses. unit : string or `~astropy.units.MagUnit` Units for the calculation, e.g., ``astropy.units.ABmag`` or ``sbpy.units.VEGAmag``. Returns ------- eff_wave : `~astropy.units.Quantity` Effective wavelengths for each ``wfb``. ci : `~astropy.units.Quantity` Color index, ``m_0 - m_1``, where 0 and 1 are element indexes for ``wfb``. """ eff_wave = [] m = np.zeros(2) * u.Unit(unit) for i in range(2): if isinstance(wfb[i], u.Quantity): if wfb[i].unit.is_equivalent(u.Hz): eff_wave.append(wfb[i].to(u.um, u.spectral())) else: eff_wave.append(wfb[i]) m[i] = self(eff_wave[i], unit=unit) elif isinstance(wfb[i], (list, tuple, SpectralElement)): w, m[i] = self.observe_bandpass(wfb[i], unit=unit) eff_wave.append(w) elif isinstance(wfb[i], str): w, pivot, m[i] = self.observe_filter_name( wfb[i], unit=unit) eff_wave.append(w) else: raise TypeError('Unsupported type for `wfb` type: {}' .format(type(wfb[i]))) ci = m[0] - m[1] return u.Quantity(eff_wave), ci
[docs]class solar_spectrum(ScienceState): """Get/set the `sbpy` default solar spectrum. To retrieve the current default: >>> from sbpy.calib import solar_spectrum >>> sun = solar_spectrum.get() To change it: >>> with solar_spectrum.set('E490_2014'): ... # E490_2014 in effect ... pass Or, you may use a string: >>> with solar_spectrum.set('E490_2014LR'): ... # E490_2014LR in effect ... pass """ _value = 'E490_2014'
[docs] @classmethod def validate(cls, value): if isinstance(value, str): return Sun.from_builtin(value) elif isinstance(value, Sun): return value else: raise TypeError( "solar_spectrum must be a string or Sun instance.")
[docs]class vega_spectrum(ScienceState): """Get/set the `sbpy` default Vega spectrum. To retrieve the current default: >>> from sbpy.calib import vega_spectrum >>> vega = vega_spectrum.get() To change it: >>> with vega_spectrum.set(Vega.from_file(filename)) # doctest: +SKIP ... # Vega from filename in effect ... pass """ _value = 'Bohlin2014'
[docs] @classmethod def validate(cls, value): if isinstance(value, str): return Vega.from_builtin(value) elif isinstance(value, Vega): return value else: raise TypeError( "vega_spectrum must be a string or Vega instance.")
class FluxdScienceState(ScienceState, ABC): """sbpy photometric calibration science states. Requires attribute: _sources. """ @classmethod def validate(cls, value): if isinstance(value, str): source = getattr(cls._sources, value) # only load the data once, save to hidden attribute ``_data`` data = source.get('data', None) if data is None: fn = get_pkg_data_filename( os.path.join('data', source['filename'])) with open(fn, 'r') as inf: phot = json.load(inf) data = {} for name, parameters in phot['data'].items(): for k, v in parameters.items(): if k == 'fluxd': k = name else: k = '{}({})'.format(name, k) data[k] = u.Quantity(v[0], v[1], subok=True) source['data'] = data return data else: return value
[docs]class solar_fluxd(FluxdScienceState): """Get/set the `sbpy` solar flux density. To set the current values: >>> from sbpy.calib import solar_fluxd >>> import sbpy.units as sbu >>> with solar_fluxd.set({'V': -26.76 * sbu.VEGAmag}): ... pass The units must be flux density per unit wavelength or per unit frequency. To retrieve the current values: >>> import sbpy.units as sbu >>> with solar_fluxd.set({'V': -26.76 * sbu.VEGAmag}): ... S = solar_fluxd.get() ... print(S['V']) -26.76 mag(VEGA) Multiple values are allowed: >>> import astropy.units as u >>> solar_fluxd.set({ ... 'PS1_g': -26.54 * u.ABmag, ... 'PS1_r': -26.93 * u.ABmag, ... 'PS1_i': -27.05 * u.ABmag ... }) # doctest: +IGNORE_OUTPUT When wavelength is required, specify ``bandpass(lambda eff)`` for effective wavelength and/or ``bandpass(lambda pivot)`` for pivot wavelength: >>> import sbpy.units as sbu >>> solar_fluxd.set({ ... 'V': -26.76 * sbu.VEGAmag, ... 'V(lambda eff)': 548 * u.nm, ... 'V(lambda pivot)': 551 * u.nm ... }) # doctest: +IGNORE_OUTPUT """ _value = 'Willmer2018' _sources = solar_sources.SolarPhotometry
[docs]class vega_fluxd(FluxdScienceState): """Get/set the `sbpy` Vega flux density. To set the current values: >>> from sbpy.calib import vega_fluxd >>> import astropy.units as u >>> with vega_fluxd.set({'V': 3674 * u.Jy}): ... pass The units must be flux density per unit wavelength or per unit frequency. To retrieve the current values: >>> import astropy.units as u >>> with vega_fluxd.set({'V': 3674 * u.Jy}): ... S = vega_fluxd.get() ... print(S['V']) 3674.0 Jy Multiple values are allowed: >>> import astropy.units as u >>> vega_fluxd.set({ ... 'PS1_g': 4026 * u.Jy, ... 'PS1_r': 3252 * u.Jy, ... 'PS1_i': 2656 * u.Jy ... }) # doctest: +IGNORE_OUTPUT When wavelength is required, specify ``bandpass(lambda eff)`` for effective wavelength and/or ``bandpass(lambda pivot)`` for pivot wavelength: >>> import astropy.units as u >>> vega_fluxd.set({ ... 'V': 3674 * u.Jy, ... 'V(lambda eff)': 548 * u.nm, ... 'V(lambda pivot)': 551 * u.nm ... }) # doctest: +IGNORE_OUTPUT """ _value = 'Willmer2018' _sources = vega_sources.VegaPhotometry
[docs]class Sun(SpectralStandard): """Solar spectrum. Most functionality requires `synphot`. The exception is retrieval of flux densities by filter name via the ``solar_fluxd`` context manager. Parameters ---------- wave : `~astropy.units.Quantity` Spectral wavelengths. fluxd : `~astropy.units.Quantity` Solar spectral flux density at 1 au. description : string, optional Brief description of the source spectrum. bibcode : string, optional Bibliography code for `sbpy.bib.register`. meta : dict, optional Any additional meta data, passed on to `~synphot.SourceSpectrum`. Attributes ---------- wave - Wavelengths of the source spectrum. fluxd - Source spectrum. description - Brief description of the source spectrum. meta - Meta data. Examples -------- Get the default solar spectrum: >>> sun = Sun.from_default() Create solar standard from `synphot.SourceSpectrum`: >>> import astropy.constants as const >>> import synphot >>> source = (synphot.SourceSpectrum(synphot.BlackBody1D, temperature=5770) ... * (3.14159 * const.R_sun**2 / const.au**2).decompose()) >>> sun = Sun(source, description='5770 K blackbody Sun') Create solar standard from a file: >>> sun = Sun.from_file('filename') # doctest: +SKIP Create solar standard from an array: >>> import numpy as np >>> import astropy.units as u >>> import astropy.constants as const >>> from astropy.modeling.blackbody import blackbody_lambda >>> wave = np.logspace(-1, 2, 300) * u.um >>> fluxd = (blackbody_lambda(wave, 5770 * u.K) * np.pi * u.sr ... * (const.R_sun**2 / const.au**2).decompose()) >>> sun = Sun.from_array(wave, fluxd, description='5770 K blackbody Sun') Interpolate to / evaluate at 0.62 μm: >>> print(sun(0.62 * u.um)) # doctest: +FLOAT_CMP 1611.7057927532567 W / (m2 um) Observe as through a spectrometer: >>> import numpy as np >>> import astropy.units as u >>> sun = Sun.from_default() >>> wave = np.linspace(1, 2.5) * u.um >>> fluxd = sun.observe(wave) # doctest: +IGNORE_OUTPUT Observe as through a filter, integrating with the bandpass transmission: >>> from sbpy.photometry import bandpass >>> sun = Sun.from_default() >>> v = bandpass('Johnson V') >>> print(sun.observe(v)) # doctest: +FLOAT_CMP 1839.93273227 W / (m2 um) Observe through a filter, using `sbpy`'s filter calibration system: >>> from sbpy.calib import solar_fluxd >>> import sbpy.units as sbu >>> solar_fluxd.set({ ... 'V': -26.76 * sbu.VEGAmag, ... 'V(lambda eff)': 548 * u.nm, ... 'V(lambda pivot)': 551 * u.nm ... }) # doctest: +IGNORE_OUTPUT >>> sun = Sun.from_default() >>> print(sun.observe('V')) -26.76 mag(VEGA) """ _sources = solar_sources.SolarSpectra _spectrum_state = solar_spectrum _fluxd_state = solar_fluxd
[docs]class Vega(SpectralStandard): """Vega spectrum. Parameters ---------- wave : `~astropy.units.Quantity` Spectral wavelengths. fluxd : `~astropy.units.Quantity` Spectral flux density. description : string, optional Brief description of the source spectrum. bibcode : string, optional Bibliography code for `sbpy.bib.register`. meta : dict, optional Any additional meta data, passed on to `~synphot.SourceSpectrum`. Attributes ---------- wave - Wavelengths of the source spectrum. fluxd - Source spectrum. description - Brief description of the source spectrum. meta - Meta data. Examples -------- Get the default Vega spectrum: >>> vega = Vega.from_default() # doctest: +IGNORE_OUTPUT Create Vega from a file: >>> vega = Vega.from_file('filename') # doctest: +SKIP Evaluate Vega at 1 μm (interpolation): >>> print(vega(1 * u.um)) # doctest: +FLOAT_CMP 6.326492514857613e-09 W / (m2 um) Observe Vega through as if through a spectrometer: >>> import numpy as np >>> wave = np.linspace(0.4, 0.6) * u.um >>> spec = vega.observe(wave) Observe Vega through a filter: >>> from sbpy.photometry import bandpass >>> V = bandpass('Johnson V') >>> fluxd = vega.observe(V) User provided calibration: >>> from sbpy.calib import vega_fluxd >>> vega = Vega(None) >>> with vega_fluxd.set({'V': 3674 * u.Jy, ... 'V(lambda pivot)': 5511 * u.AA}): ... print(vega.observe('V', unit='Jy')) 3674.0 Jy """ _sources = vega_sources.VegaSpectra _spectrum_state = vega_spectrum _fluxd_state = vega_fluxd