Source code for species.read.read_filter

"""
Module with reading functionalities for filter profiles.
"""

import os
import warnings

from configparser import ConfigParser
from typing import Union, Tuple

import h5py
import numpy as np

from typeguard import typechecked
from scipy.interpolate import interp1d, InterpolatedUnivariateSpline

from species.data.filter_data.filter_data import add_filter_profile
from species.data.spec_data.spec_vega import add_vega


[docs] class ReadFilter: """ Class for reading a filter profile from the database. """ @typechecked def __init__(self, filter_name: str) -> None: """ Parameters ---------- filter_name : str Filter name as stored in the database. Filter names from the SVO Filter Profile Service will be automatically downloaded, stored in the database, and read from the database. Returns ------- NoneType None """ self.filter_name = filter_name config_file = os.path.join(os.getcwd(), "species_config.ini") config = ConfigParser() config.read(config_file) self.database = config["species"]["database"] self.data_folder = config["species"]["data_folder"] with h5py.File(self.database, "r") as hdf5_file: # Check if the filter is found in 'r' mode # because the 'a' mode is not possible when # using multiprocessing if f"filters/{self.filter_name}" in hdf5_file: filter_found = True else: filter_found = False if not filter_found: with h5py.File(self.database, "a") as hdf5_file: add_filter_profile(self.data_folder, hdf5_file, self.filter_name)
[docs] @typechecked def get_filter(self) -> np.ndarray: """ Select a filter profile from the database. Returns ------- np.ndarray Array with the wavelengths and filter transmission. The array has 2 dimensions with the shape (n_wavelengths, 2). """ with h5py.File(self.database, "r") as hdf5_file: data = np.asarray(hdf5_file[f"filters/{self.filter_name}"]) if data.shape[0] == 2 and data.shape[1] > data.shape[0]: # Required for backward compatibility data = np.transpose(data) return data
[docs] @typechecked def interpolate_filter(self) -> interp1d: """ Interpolate a filter profile with the `interp1d <https:// docs.scipy.org/doc/scipy/reference/generated/ scipy.interpolate.interp1d.html>`_ function from ``scipy.interpolate`` and linear kind of interpolation. Returns ------- scipy.interpolate.interp1d Linearly interpolated filter profile. """ data = self.get_filter() return interp1d( data[:, 0], data[:, 1], kind="linear", bounds_error=False, fill_value=float("nan"), )
[docs] @typechecked def wavelength_range( self, ) -> Tuple[Union[np.float32, np.float64], Union[np.float32, np.float64]]: """ Extract the wavelength range of the filter profile. Returns ------- float Minimum wavelength (:math:`\\mu\\mathrm{m}`). float Maximum wavelength (:math:`\\mu\\mathrm{m}`). """ data = self.get_filter() return data[0, 0], data[-1, 0]
[docs] @typechecked def mean_wavelength(self) -> Union[np.float32, np.float64]: """ Calculate the weighted mean wavelength of the filter profile. Returns ------- float Mean wavelength (:math:`\\mu\\mathrm{m}`). """ data = self.get_filter() return np.trapz(data[:, 0] * data[:, 1], x=data[:, 0]) / np.trapz( data[:, 1], x=data[:, 0] )
[docs] @typechecked def effective_wavelength(self) -> Union[np.float32, np.float64]: """ Calculate the effective wavelength of the filter profile. The effective wavelength is calculated as the weighted average based on the filter profile and the spectrum of Vega. Returns ------- float Effective wavelength (:math:`\\mu\\mathrm{m}`). """ filter_profile = self.get_filter() with h5py.File(self.database, "a") as hdf5_file: if "spectra/calibration/vega" not in hdf5_file: add_vega(self.data_folder, hdf5_file) vega_spec = np.array(hdf5_file["spectra/calibration/vega"]) flux_interp = interp1d( vega_spec[0,], vega_spec[1,], bounds_error=False, fill_value="extrapolate", ) flux_filter = flux_interp(filter_profile[:, 0]) return np.trapz( filter_profile[:, 0] * filter_profile[:, 1] * flux_filter, x=filter_profile[:, 0], ) / np.trapz(filter_profile[:, 1] * flux_filter, x=filter_profile[:, 0])
[docs] @typechecked def filter_fwhm(self) -> float: """ Calculate the full width at half maximum (FWHM) of the filter profile. Returns ------- float Full width at half maximum (:math:`\\mu\\mathrm{m}`). """ data = self.get_filter() spline = InterpolatedUnivariateSpline( data[:, 0], data[:, 1] - np.max(data[:, 1]) / 2.0 ) root = spline.roots() diff = root - self.mean_wavelength() root1 = np.amax(diff[diff < 0.0]) root2 = np.amin(diff[diff > 0.0]) return root2 - root1
[docs] @typechecked def effective_width(self) -> Union[np.float32, np.float64]: """ Calculate the effective width of the filter profile. The effective width is equivalent to the horizontal size of a rectangle with height equal to the maximum transmission and with the same area as the one covered by the filter profile. Returns ------- float Effective width (:math:`\\mu\\mathrm{m}`). """ data = self.get_filter() return np.trapz(data[:, 1], x=data[:, 0]) / np.amax(data[:, 1])
[docs] @typechecked def detector_type(self) -> str: """ Return the detector type. Returns ------- str Detector type ('energy' or 'photon'). """ with h5py.File(self.database, "r") as hdf5_file: dset = hdf5_file[f"filters/{self.filter_name}"] if "det_type" in dset.attrs: det_type = dset.attrs["det_type"] else: warnings.warn( f"Detector type not found for {self.filter_name}. " "The database was probably created before the " "detector type was introduced in species v0.3.1. " "Assuming a photon-counting detector." ) det_type = "photon" return det_type