import yaml
import os
import os.path
import math
import numpy as np
import pickle
import logging
from typing import Dict, Tuple, List, Union, Any
Show_Configuration = Dict[ str, Dict[str, Union[float,int,str,List[float],List[int]]] ]
#------------------------------------------------------------------------------
# read_configuration
#------------------------------------------------------------------------------
[docs]def read_configuration( anycase_instrumentname : str) ->Show_Configuration:
instrumentname = anycase_instrumentname.lower() # Get the instrument name in lower case
fullname = __file__ # Get the filename of this module
[basedir, spare] = os.path.split(fullname) # split out the list_of_files
site_specific_configuration = {} # load in site specific settings
fullname = os.path.normpath(basedir + os.sep + 'yaml_files' + os.sep + instrumentname + '_site_settings.yaml')
try:
f = open(fullname, 'rt')
site_specific_configuration = yaml.load(f)
f.close()
site_specific_configuration['files']['site_specific_filename']=fullname
except Exception as inst:
print(type(inst)) # the exception instance
print(inst.args) # arguments stored in .args
print(inst) # __str__ allows args to be printed directly,
raise RuntimeError('Error reading site specific configuration file for instrument %s from file %s'%(anycase_instrumentname,fullname))
basename = site_specific_configuration['files'].get('cdb_base_directory')
if (basename is None):
fullname = os.path.normpath(basedir + os.sep + 'yaml_files' + os.sep + instrumentname + '_default.yaml')
else:
fullname = os.path.normpath(basename + os.sep + instrumentname + '.yaml')
configuration = {}
try:
f = open(fullname, 'rt')
configuration = yaml.load(f)
f.close()
except Exception as inst:
print(type(inst)) # the exception instance
print(inst.args) # arguments stored in .args
print(inst) # __str__ allows args to be printed directly,
#raise RuntimeError('Error reading instrument configuration for instrument %s from file %s'%(anycase_instrumentname,fullname))
for k,v in site_specific_configuration.items():
configuration[k] = v
return configuration
#------------------------------------------------------------------------------
# decode_version
#------------------------------------------------------------------------------
[docs]def decode_version( versionstr: str):
fields = versionstr.split('.')
version = int(fields[0]), int(fields[1]), int(fields[2])
return version
#------------------------------------------------------------------------------
# encode_version
#------------------------------------------------------------------------------
[docs]def encode_version( version: Tuple[int,int,int])->str:
versionstr = "%02d.%02d.%02d"%(version[0], version[1], version[2])
return versionstr
#------------------------------------------------------------------------------
# cdb_directory
#------------------------------------------------------------------------------
[docs]def cdb_directory( config :Dict[str,Any]) ->str :
"""
Returns the location of the instrument Calibration Database Directory
:param anycase_instrumentname:
:return:
"""
fullname = config["files"]["cdb_base_directory"]
return (os.path.normpath(fullname)) + os.sep
#------------------------------------------------------------------------------
# read_pickle_h20xsection
#------------------------------------------------------------------------------
[docs]def read_pickle_h20xsection() -> Tuple[np.ndarray, np.ndarray]:
fullname = __file__ # Get the filename of this module
[basedir, spare] = os.path.split(fullname) # split out the list_of_files
fullname = basedir + os.sep + 'ancillary_data_files' + os.sep + 'pickle_h2o_xsection.dat'
try:
f = open(fullname, 'rb')
w,absxs = pickle.load(f)
f.close()
except Exception as inst:
raise RuntimeError('Error reading water cross-section data from file configuration file from file %s'%(fullname,))
return w,absxs
#------------------------------------------------------------------------------
# Parameters
#------------------------------------------------------------------------------
[docs]class Parameters:
"""
The analysis of SHOW involves transformations from interferogram
space to frequency space. This class provides an interface that
converts between the two spaces using known and/or fitted design
parameters.
The interferogram equation is given by:
I(x) = 0.5*Integral( B(k)(1 +cos(2.pi.k.x)dk ) (eqn 1)
where k is the heterodyne frequency:
k = 4M( v - vo)tan(ThetaL) (eqn 2)
or (v-v0) = k/(4Mtan(thetaL))
M = Magnification (typically 4 to 5)
v = wavenumber cm-1
vo = Littrow wavenumber cm-1
ThetaL = Littrow angle of grating (typically ~28.5 degrees).
The measured interferogram is imaged onto a detector of width 'L' cms
and the discrete harmonics of equation(1) that are retrieved by a
discrete Fourier Transform are given by:
k.L = n n = 1,2,3,4,5,6,7,8.... (eqn 3)
and by substitution from eqn (2)
n
(vn - v0) = ------------------ (eqn 4)
4.M.L.tan(thetaL)
where vn is the nth discrete wavenumber.
"""
#----------------------------------------------------------------
# Constructor::Parameters
#----------------------------------------------------------------
def __init__(self, instrumentname : str = "er2_2017") -> None :
"""
:param instrumentname:
The name of the instrument configuration. This will be the name part of a yaml file stored in the config
subdirectory. The following values were supported as of 2017-01-10
- er2_2017
- timmins_2014
- uofs_dec2016
The following fields are defined
- config, The entire dictionary read in from the YAML file.
- instrumentname, The name of the instrument (and the YAML file).
- littrow_angle The grating littrow angle in radians.
- ggd, grating groove density, lines per cm
- magnification, magnification of the back end optics
- num_spectralpixels, the number of detector pixels in the spectral direction
- spectralpixel_width_cms, the width of each pixel in the spectral direction in cms.
- fft_window_start, the (inclusive) lower index in the spectral dimension for the spectral sub-window given to the FFT algorithm
- fft_window_end, the (inclusive) upper index in the spectral dimension for the spectral sub-window given to the FFT algorithm
- height_window_start the (inclusive) lower index in the height dimension used to generate the level 1 sub-window. We only need the useful area.
- height_window_end the (inclusive) upper index in the height dimension used to generate the level 1 sub-window. We only need the useful area.
"""
self.config = read_configuration( instrumentname ) # type: Dict[str, Any]
# self.config = read_subwindow_configuration( self.config)
self.instrumentname = instrumentname # type: str
config = self.config['shsdesign'] #
self.littrow_angle = math.radians(config['littrow_angle']) # type: float
self.ggd = config['grating_groove_density_percm'] # type: float
self.magnification = 1.0/self.config['cdb']['optics_magnification'] # type: float
self.num_spectralpixels = int( config['num_spectralpixels']) # type: int
self.spectralpixel_width_cms = config['spectralpixel_width_cms'] # type: float
self.fft_window_start = int( self.config['iFOV']['Pixel_column_at_left_edge_of_grating_image']) # type: int
self.fft_window_end = int( self.config['iFOV']['Pixel_column_at_right_edge_of_grating_image']) # type: int
self.height_window_start = int( self.config['iFOV']['Pixel_row_at_bottom_of_atmospheric_image']) # type: int
self.height_window_end = int( self.config['iFOV']['Pixel_row_at_top_edge_of_atmospheric_image']) # type: int
self.ccd_width_cms = (self.fft_window_end - self.fft_window_start +1) * self.spectralpixel_width_cms # type: float
# These parameters change as the instrument temperature changes
#self.littrow_wavelen_nm = config['littrow_wavelen_nm'] # Littrow wavelength nm, e.g. 1363.44018518519
#self.littrow_wavenumber = 1.0E7 / self.littrow_wavelen_nm # Littrow wavenumber cm-1,
#------------------------------------------------------------------------------
# Parameters::configuration_section
#------------------------------------------------------------------------------
[docs] def section( self, section: str) -> Dict :
"""
Returns a specific section of the configuration parameters
- shsdesign
- level0
The code will throw an exception if the requested section does not exist.
:param section:
Name of the configuration section to return
:return:
The configuration ssection
"""
return self.config[section]
# ------------------------------------------------------------------------------
# Parameters::active_spectral_length_cms
# ------------------------------------------------------------------------------
[docs] def active_length_cms(self):
return (self.fft_window_end - self.fft_window_start) * self.spectralpixel_width_cms
#------------------------------------------------------------------------------
# fft_step_size_wavenumber
#------------------------------------------------------------------------------
[docs] def fft_step_size_wavenumber( self ) -> float:
"""
Return the theoretical step size in wavenumbers of one
harmonic from the discrete FFT.
:return: The spacing of one FFT harmonic in wavenumbers (cm-1)
:rtype: float
"""
M = self.magnification
L = self.active_length_cms()
tantheta = math.tan( self.littrow_angle)
dsigma = 1.0/( 4.0*M*L*tantheta)
return dsigma
#------------------------------------------------------------------------------
# gauss_pixel_sigma_to_wavenumber_sigma
#------------------------------------------------------------------------------
[docs] def gauss_pixel_sigma_to_wavenumber_sigma(self, pixelsigma : (float, np.ndarray)) -> (float, np.ndarray) :
"""
Calculate the standard deviation width of a Gaussian in
FFT wavenumber space given the sd width of the corresponding gaussian
in CCD pixel space. This is useful if we want to multiply in
detector space but convolve theoretical spectra in wavenumber space.
:param pixelsigma
:return: The standard deviation width of the gaussian in wavenumbers cm-1.
:rtype: float or numpy.ndarray
"""
M = self.magnification
sigmacms = pixelsigma * self.active_length_cms() # Get the width of the Gaussian in centimeters on the detector
ksigma = 1.0/(2*math.pi*sigmacms) # Get the width of the Gaussian in 'k' space
wavenumbersigma = ksigma/(4 * M * math.tan( self.littrow_angle)) # Get the width of the Gaussian in wavenumbers in FFT space
return wavenumbersigma
#------------------------------------------------------------------------------
# gauss_wavenumber_sigma_to_pixel_sigma
#------------------------------------------------------------------------------
[docs] def gauss_wavenumber_sigma_to_pixel_sigma(self, wavenumbersigma:(float, np.ndarray)) -> (float, np.ndarray):
"""
Calculate the standard deviation width of a Gaussian in
CCD pixel space given the sd width of the corresponding gaussian
in FFT wavenumber space. This is useful if we want to multiply in
wavenumber space but convolve in CCD pixel space.
:return: The sd width of the gaussian in CCD pixels .
:rtype: float or numpy.ndarray
"""
M = self.magnification # Get the magnification
ksigma = wavenumbersigma*4*M*math.tan(self.littrow_angle) # Get the width of the Gaussian in 'k' space
sigmacms = 1.0/(2*math.pi*ksigma) # Get the width of the Gaussian in CCD space in centimeters
pixelsigma = sigmacms / self.spectralpixel_width_cms # Get the width of the Gaussian in CCD space in pixels
return pixelsigma
#------------------------------------------------------------------------------
# Gauss_PixelSigmaToWavelenSigma
#------------------------------------------------------------------------------
[docs] def gauss_pixel_sigma_to_wavelen_sigma(self, pixelsigma :(float, np.ndarray)) -> (float, np.ndarray):
"""
Calculate the standard deviation width of a Gaussian in
FFT wavelength space given the sd width of the corresponding gaussian
in CCD pixel space. This is useful if we want to multiply in
detector space but convolve theoretical spectra in wavenumber space.
:return: The sd width of the gaussian in wavenumbers cm-1.
:rtype: float or numpy.ndarray
"""
wavenumbersigma = self.gauss_pixel_sigma_to_wavenumber_sigma(pixelsigma)
lam1 = 1.0E7 /self.littrow_wavenumber
lam2 = 1.0E7/(self.littrow_wavenumber + wavenumbersigma)
wavelensigma = lam1-lam2
return wavelensigma
#------------------------------------------------------------------------------
# nominal_littrow_wavenumber
#------------------------------------------------------------------------------
g_onetime = 0
[docs] def nominal_littrow_wavenumber(self, shs_temperature: float ):
if (Parameters.g_onetime < 1):
logging.warning('Parameters.nominal_littrow_wavenumber, still needs the temperature dependence added in properly')
Parameters.g_onetime += 1
littrow = self.config['cdb']['nominal_littrow_wavelen_nm']
mu0 = 1.0E7/littrow # for the moment insert a nominal estimate. Later on replace with temperature dependence
return mu0
#------------------------------------------------------------------------------
# NominalWavenumbers
#------------------------------------------------------------------------------
[docs] def nominal_wavenumbers(self, numfreq:int, shs_temperature_celsius: float)-> np.ndarray:
"""
Calculates the nominal wavenumbers of the retrieved spectrum
for the first "numfreq" discrete harmonics.
:param numfreq: The number of frequencies.
:type numfreq: int
:return: The array (numfreq) of nominal wavenumbers in cm-1
:rtype: numpy.ndarray
"""
deltasigma = self.fft_step_size_wavenumber()
freq = np.arange(0,numfreq)
nominal_wavenum = self.nominal_littrow_wavenumber(shs_temperature_celsius) - freq * deltasigma # A rough and ready wavelength calibration. No stretch and shift yet
return nominal_wavenum
#------------------------------------------------------------------------------
# NominalWavelengths
#------------------------------------------------------------------------------
[docs] def nominal_wavelengths(self, numfreq:int, shs_temperature_celsius: float)-> np.ndarray:
"""
Calculates the nominal wavelengths of the retrieved spectrum
for the first "numfreq" discrete harmonics.
:param numfreq: The number of frequencies.
:type numfreq: int
:return: The array (numfreq) of nominal wavenumbers in cm-1
:rtype: numpy.ndarray
"""
nominal_wavelen = self.nominal_wavenumbers(numfreq, shs_temperature_celsius)
nominal_wavelen = 1.0E7/nominal_wavelen
return nominal_wavelen