This notebook forms part of a series on computational optical radiometry. The notebooks can be downloaded from Github. These notebooks are constantly revised and updated, please revisit from time to time.

This notebook was last tested with Python 3.6.1. Consider moving to Python 3.6 if you are not already using it.

The date of this document and module versions used in this document are given at the end of the file.
Feedback is appreciated: neliswillers at gmail dot com.

## Overview

This notebook provides a description of the basics of radiometric spectral variables and conversions between spectral densities, generic filter functions, generic detector functions, and reading and plotting spectral data. It then moves on to how to calculate spectral integrals, calculate spectral effective values, spectral convolution, colour coordinate calculations, and spatial integrals.



In [1]:

from IPython.display import display
from IPython.display import Image
from IPython.display import HTML



## Spectral variables

Consider the mathematical sine function $s(t) = A \sin(2\pi f t)$. Here $f$ is the frequency, $t$ is time and $A$ is the amplitude. As a mathematical function all three these parameters determine the value of the function. A sine wave generator on an engineer's workbench creates a temporal sine wave, i.e., a wave shape of sinusoidal characteristics, that varies with time. The engineer can adjust the properties of the signal by changing the frequency $f$ or amplitude $A$ by turning a knob on the instrument. The engineer cannot adjust the time $t$, because time is not a free variable in the sine wave generator, time flows as part of the human reality. We can say that the sine wave depends only on $f$ and $A$, but it varies with time $t$.

Imagine a pipe with water flowing through it. One can measure the volume of water flowing past a fixed point in a given time period; the longer the time period, the larger the water volume. It makes sense to define a temporal flow rate expressed e.g., in cubic metres per second (at given time $t$). Now generalise this notion to some form of density with units of [something/s], [something/m$^3$], or [something/$\mu$m] These variables can be integrated over time/volume/spectral density to arrive at a cumulative quantity. The key observation here is that a quantity can be expressed as a density over time or volume, or in the radiometry case a quantity over a spectral band.

The spectral domain is similar to the time domain: it is not a free variable that can be adjusted, but rather provides a domain of observation against which an object's properties can be calculated or measured. These properties are said to vary with the spectral variable. Some spectral properties depend on free variables, e.g., Planck radiance depends on the temperature of the object, but varies with spectral variable.

Spectral domains (Sec 2.3.3) are commonly expressed in wavelength $\lambda$, wavenumber $\tilde{\nu}$ (the number of waves that will fit into a 1-cm length, written as one word) or frequency $f$. In the optical domain wavelength and wavenumber are more commonly used than frequency. Wavelength is commonly given in units of nanometres (nm) or micrometres ($\mu$m), wavenumber in (cm$^{-1}$) and frequency in [Hz].

Spectral variables come in two varieties: variable without density (equivalent to the sine wave) and density variables (equivalent to the water in the pipe). Variables without density include spectral transmittance of a filter, the atmosphere or emissivity. Variables with density include Planck radiation. Spectral density values have units of [something/$\mu$m], [something/cm$^{-1}$] or [something/Hz].

The total quantity of a spectral variable over a spectral range is determined by integrating the spectral variable over the spectral range, similar to the water example above.

Wavelength and wavenumber are related by $\lambda = 10^4/\tilde{\nu}$ (Sec 2.3.3). The conversion of spectral densities requires the derivative of the previous relation $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$.

Spectral variables are represented in Python as Numpy arrays (Sec D.3). These arrays are rank-one (N,) or rank-two (N,1), but always represent the spectral variable as a vector of values. The conversion between rank one and rank two is done by reshaping the vector.



In [2]:

%matplotlib inline
import numpy as np

#create a micrometre wavelength spectral range with 101 samples
#from 400 nm to 800 nm
wl1 = np.linspace(0.4, 0.8, 101) #rank (N,)
wl2 = wl1.reshape(-1,1) #convert to rank (N,1)
wl3 = wl2.reshape(1,-1) #convert to rank (1,N)
wl4 = wl2.reshape(-1,)  #convert back to rank (N,)
print('Wavelength vector as rank one, has shape {}'.format(wl1.shape))
print('Wavelength vector as rank two, has shape {}'.format(wl2.shape))
print('Wavelength vector as rank two, has shape {}'.format(wl3.shape))
print('Wavelength vector as rank one, has shape {}'.format(wl4.shape))




Wavelength vector as rank one, has shape (101,)
Wavelength vector as rank two, has shape (101, 1)
Wavelength vector as rank two, has shape (1, 101)
Wavelength vector as rank one, has shape (101,)



It appears that in Numpy a rank-one array (N,) behaves like a row vector (1,N) in some instances. In pyradi it is more convenient and natural to consider spectral variables as column vectors (N,1) to this rehaping from (N,) to (N,1) takes place quite frequently. Why consider spectral variables as (N,1) arrays? When saving spectral data to file, it seems more natural (and easier to view/edit) if the spectral range varies along the rows in the file: each spectral value is on a new line in the file. This approach also lends itself to saving more than one spectral variable in a file where the first column represents the spectral domain and all remaining columns represent spectral variables. A case in point is the Modtran tape7 file. So as a general rule: pyradi always expect spectral variables of one of the these shapes:



In [3]:

wl1 = np.linspace(0.4, 0.8, 101) #rank (N,)
wl2 = wl1.reshape(-1,1) #convert to rank (N,1)



The pyradi library provides several conversion functions to convert between different spectral domains and also to convert between spectral densities. The units must consistently be in micrometres ($\mu$m), wavenumber in (cm$^{-1}$) or frequency in [Hz] - any other units will provide incorrect results. These functions are documented here pyradi.ryutils.convertSpectralDomain and pyradi.ryutils.convertSpectralDensity.



In [4]:

wl = np.linspace(1, 14, 101).reshape(-1,1) # wavelength
wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber
f = ryutils.convertSpectralDomain(wn, type='nf') # frequency

p = ryplot.Plotter(1,1,3,figsize=(12,4));
p.plot(1,wl,wn,'','Wavelength $\mu$m', 'Wavenumber cm$^{-1}$');
p.plot(2,wl,f,'','Wavelength $\mu$m', 'Frequency Hz');
p.plot(3,wn,f,'','Wavenumber cm$^{-1}$', 'Frequency Hz',maxNX=4);






Converting spectral density values are similar to converting between spectral domains. In this example the pyradi.ryplanck.planck is used to calculate the spectral density of Planck radiance (Sec 3.1, Sec D.4.1) for a temperature of 1000 K as a wavelength spectral density. This value is then converted to a wavenumber spectral density. Finally, the Planck radiance is calculated in the wavenumber domain, which should give the same result as the converted spectral radiance.



In [5]:

wl = np.linspace(1, 14, 101)#.reshape(-1,1) # wavelength
wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber

radiancewl = ryplanck.planck(wl,1000.0, 'el') / np.pi # W/(m2.sr.um)
radiancewn2 = ryplanck.planck(wn,1000, 'en') / np.pi

p = ryplot.Plotter(2,1,3,figsize=(12,4))
p.plot(1,wl, radiancewl,'Planck radiance','Wavelength $\mu$m',
'Radiance W/(m$^2$.sr.$\mu$m)',maxNX=4)
p.plot(2,wn, radiancewn1,'Planck radiance','Wavenumber cm$^{-1}$',
'Radiance W/(m$^2$.sr.cm$^{-1}$)',maxNX=4)
p.plot(3,wn, radiancewn2,'Planck radiance','Wavenumber cm$^{-1}$',
'Radiance W/(m$^2$.sr.cm$^{-1}$)',maxNX=4);






See the section Performing spectral integrals below on performing integrals on spectral data.

## Interpolation of Spectral Variables

SciPy has a number of interpolation options in scipy.interpolate. The one-dimensional interpolation scheme enables the definition of the lookup table as a function by using scipy.interpolate.interp1d. The user can define the interpolation scheme to be used - the default is kind='linear', to do cubic interpolation use kind='cubic'. In the example below observe the difference between the two schemes.

The example below is is an ill-posed problem: normally one would have a much denser input data set. The purpose with using this poor data set it to illustrate the dangers of interpolation. Do not use a potentially unstable interpolation scheme such as cubic interpolation on a poor data set such as the one shown below

The input data set has only four samples. The ordinate (y) data has a slight oscillation. The abscissa (x) data are defined in the wavelength space. Equivalently, the abscissa (x) data are also defined in the wavenumber space, because $\tilde{\nu}=10^4/\lambda$. This is the only information available, there is no requirement on the values between the stated abscissa.

Important: by selecting the interpolation scheme (linear, cubic, etc.), you are implicitly making a modelling assumption. Linear interpolation in the wavelength space means a non-linear interpolation in wavenumber space and vice versa (because wavelength and wavenumber are nonlinearly related). Stated differently: if your input data set was compiled to be interpolated linearly in wavelength space (blue lines), then linear interpolation in wavenumber space will give the wrong results (red lines).

Also, by selecting between the linear and cubic (or any other scheme) you are also making a modelling assumption. Compare the difference between the results for the linear interpolation (red/blue) and the cubic interpolation (green/cyan) lines.

It is evident in the cubic interpolation graphs (green/cyan) that cubic/high-order interpolation schemes do not take kindly to large variations in abscissa increments. For example, the large increment from 10 to 20 in the wavelength space (green line) resulted in a huge variation in the spline curve. Likewise, the large wavenumber increment from 5000 to 10000 (or 1 to 2 $\mu$m) (cyan line) causes a large variation.

From this experiment we can draw the following conclusions:

1. Clearly state the space (wavelength/wavenumber) where the input data set applies. Do interpolation only in the prescribed space.

2. Ensure that the input data set samples at similar/even density throughput the abscissa domain. If no input data is available in some subranges in the abscissa domain, construct synthetic data by some modelling assumptions, to ensure even density. Such an even spread in abscissa values will limit wide swings in the high-order interpolation schemes.

Thanks to Riana Willers for pointing out the subtle matter of interpolation in wavelength or wavenumber space.



In [6]:

from scipy import interpolate
import numpy as np
%matplotlib inline

# create input data in wavelength and wavenumber domain
maxx = 24
wl = np.asarray([1, 2, 5, 10, maxx])
wn = 1e4 / wl # wl and wn values precisely align
fn = np.asarray([2,3,1, 2, -1])

# create the plotting samples - much tighter than input data
wli = np.linspace(1,maxx,400)
wni = 1e4 / wli

#create the lookup functions
fwll = interpolate.interp1d(wl,fn)
fwlc = interpolate.interp1d(wl,fn, kind='cubic')
fwnl = interpolate.interp1d(wn,fn)
fwnc = interpolate.interp1d(wn,fn, kind='cubic')

p = ryplot.Plotter(1,2,2,figsize=(12,8))
p.plot(1,wli,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(1,wli,fwnl(wni),label=['wn scale interpolation - linear'],plotCol=['r'], maxNX=5)
p.resetPlotCol()
p.plot(2,wli,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(2,wli,fwlc(wli),label=['wl scale interpolation - cubic'],plotCol=['g'])
p.plot(2,wli,fwnc(wni),label=['wn scale interpolation - cubic'],plotCol=['c'], maxNX=5)
p.resetPlotCol()
p.plot(3,wni,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(3,wni,fwnl(wni),label=['wn scale interpolation - linear'],plotCol=['r'], maxNX=5)
p.resetPlotCol()
p.plot(4,wni,fwll(wli),label=['wl scale interpolation - linear'],plotCol=['b'])
p.plot(4,wni,fwlc(wli),label=['wl scale interpolation - cubic'],plotCol=['g'])
p.plot(4,wni,fwnc(wni),label=['wn scale interpolation - cubic'],plotCol=['c'], maxNX=5);






## Generic spectral filter response

In most real-world system analyses, filter spectral responses will be measured and used in tabular formats. However, in some instances it is convenient to create spectral filter shapes according to a parameter set. The pyradi library has a function pyradi.ryutils.sfilter to assist in calculating a variety of spectral filter shapes. This function has the form (Sec D.4.2):

$$\tau_\lambda =\tau_s + \tau_p \exp\left[ -\left(\frac{2(\lambda-\lambda_c)}{\Delta\lambda}\right)^s\right].$$

This function takes a spectral variable $\lambda$, the centre wavelength $\lambda_c$, the spectral width $\Delta\lambda$, the sharpness of cutoff $s$ (if $s=2$, the curve has a Gaussian shape, and if $s=\infty$, the curve has a square top-hat shape), transmittance in the pass band $\tau_p$ and stop band $\tau_s$ and the type of filter. The use of the function is best demonstrated by graphical means.



In [7]:

wavelength = np.linspace(0.1, 1.3, 350).reshape(-1, 1)

#demonstrate the different filter types
width = 0.4
center = 0.9
filterExp=[2, 2,  4, 6,  8, 12, 500] # shaprness
filterPass=[0.4, 0.5,  0.6, 0.7,  0.8, 0.9, 0.99] # tau pass
filterSupp = np.asarray(filterPass) * 0.1 # tau stop
filterType=['bandpass', 'lowpass', 'highpass', 'bandpass', 'lowpass',
'highpass', 'bandpass']
filterTxt = [r's={0}, $\tau_p$={2}, {1}'.format(s,y,z) for s,y,
z in zip(filterExp, filterType, filterPass) ]
filters = ryutils.sfilter(wavelength,center, width, filterExp[0],
filterPass[0],  filterSupp[0], filterType[0])

for i,exponent in enumerate(filterExp[1:]):
tau=ryutils.sfilter(wavelength,center, width, filterExp[i+1],filterPass[i+1],
filterSupp[i+1], filterType[i+1])
filters =  np.hstack((filters,tau))

#plot sample filters
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4));
smpleplt.plot(1, wavelength, filters,
r"Optical filter for $\lambda_c$={0}, $\Delta\lambda$={1}".format(center,width),
r'Wavelength [$\mu$m]', r'Transmittance', \
['r', 'g', 'y','g', 'b', 'm','c'],label=filterTxt);







In [8]:

#all passband, different shapes
width = 0.5
center = 0.7
filterExp=[2,  4, 6,  8, 12, 500] # sharpness

filterTxt = ['s={0}'.format(s) for s in filterExp ]
filters = ryutils.sfilter(wavelength,center, width, filterExp[0], 0.8,  0.1)
for exponent in filterExp[1:]:
filters =  np.hstack((filters, ryutils.sfilter(wavelength,center, width,
exponent, 0.8,  0.1)))
smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4));
smpleplt.plot(1, wavelength, filters,
r"Optical filter for $\lambda_c$=0.7,$\Delta\lambda$=0.5,$\tau_{s}$=0.1, $\tau_{p}$=0.8",
r'Wavelength [$\mu$m]', r'Transmittance', \
['r', 'g', 'y','g', 'b', 'm'],label=filterTxt);






## Detector spectral response

Detector spectral responsivity values are also normally measured or taken from detector data sheets. In the absence of such information the shape can be calculated from a detailed bulk detector model, such as provided in pyradi.rydetector, but sometimes a simplified shape is required. Such a simple photon detector responsivity can be calculated with the pyradi.ryutils.responsivity function. The function takes the form (Sec D.4.3):

$${\cal R}_\lambda= k\left[ \left(\frac{\lambda}{\lambda_c} \right)^a - \left(\frac{\lambda}{\lambda_c} \right)^n \right]\;\;\;{\rm for}\;\;\; {\cal R}_\lambda>0, \;\;\; 0 \;\;\;{\rm otherwise}.$$

This function takes a wavelength spectral variable $\lambda$, the peak wavelength value $\lambda_c$, the sharpness of cuton (lower wavelength range $a$, the sharpness of cutoff (longer wavelength range) $n$, and an amplitude scaling factor $k$. Best results are obtained for $0\leq a \leq 5$ and $5\leq n \leq 50$. The use of the function is best demonstrated by graphical means.



In [9]:

lwavepeak = 1.2
params = [(0.5, 5), (1, 10), (1, 20), (1, 30), (1, 1000), (2, 20)] #cut-on and cut-off
parameterTxt = ['a={0}, n={1}'.format(s[0], s[1]) for s in params ]
responsivities = ryutils.responsivity(wavelength,lwavepeak, params[0][0], params[0][1], 1.0)
for param in params[1:]:
responsivities =  np.hstack((responsivities, ryutils.responsivity(wavelength,lwavepeak, param[0], param[1], 1.0)))

smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4))
smpleplt.plot(1, wavelength, responsivities, "Detector Responsivity for $\lambda_c$=1.2 $\mu$m, k=1", r'Wavelength [$\mu$m]',\
r'Responsivity', \
['r-', 'g-', 'y-','g--', 'b-', 'm-'],label=parameterTxt);







In [10]:

# here we simply multiply the responsivity and spectral filter spectral curves.
# this is a silly example, but demonstrates the spectral integral.
filtreps = responsivities * filters
parameterTxt = [str(s)+' & '+str(f) for (s, f) in zip(params, filterExp) ]

smpleplt = ryplot.Plotter(1, 1, 1, figsize=(10, 4))
smpleplt.plot(1, wavelength, filtreps, "Filtered Detector Responsivity",
r'Wavelength $\mu$m',\
r'Responsivity', \
['r-', 'g-', 'y-','g--', 'b-', 'm-'],label=parameterTxt,legendAlpha=0.5);






## Plotting of spectral variables

The title appears trivial, but there is a hidden subtlety when working with and plotting radiometric variables. In the section above spectral variables in pyradi, the data was created in the wavelength domain (using the linspace function), but plotted in both wavelength and wavenumber domains. Note that the data was not resampled to wavenumber domain before plotting: the wavenumber plots are not at constant intervals. The wavelength domain was merely recalculated to wavenumber domain. This is illustrated in the graph below.



In [11]:

wlc = np.linspace(0.5, 8.5, 17).reshape(-1,1) # wavelength
wnc = ryutils.convertSpectralDomain(wlc, type='ln') # wavenumber
p = ryplot.Plotter(2,1,3,figsize=(18,2))
for i in range(wlc.shape[0]):
p.plot(1,np.array([wlc[i],wlc[i]]), np.array([0,1.]), 'Constant wavelength intervals',
'Wavelength $\mu$m', plotCol=['r'], markers=['o'], maxNY=1)
p.plot(2,np.array([wnc[i],wnc[i]]), np.array([0,1.]), 'Non-constant wavenumber intervals',
'Wavenumber cm$^{-1}$',plotCol=['r'], markers=['o'], maxNX=5, maxNY=1, pltaxis=[0, 10000, 0, 1]);






The data and the plot are equally 'valid' or true in either domain, because the domain merely samples the data at different spectral intercepts. So we can load and process data in either domain and plot in either domain, without losing any truth value. Plot in whatever domain illustrates your argument the best.

You should however take care when dealing with spectral densities: it is possible and mathematically allowable to plot a wavenumber density with units of W/(m$^2$.sr.cm$^{-1}$) against a wavelength domain with units of $\mu$m, but it does not make much practical sense to do so.

Consider the intervals in both cases, here calculated with numpy.diff. Note that the intervals are constant and positive for the wavelength scale, but negative and non-constant for the wavenumber scale as converted from the wavelength scale. Recall from Section Spectral-variables above, that $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$. In other words the spectral intervals change in the opposite sense: the wavelength domain scale increases (positive intervals), but when converted to wavenumber domain the scale decreases (negative intervals).



In [12]:

p = ryplot.Plotter(3,1,3,figsize=(18,2))
p.plot(1,(wlc[1:]+wlc[:-1])/2, np.diff(wlc,axis=0), 'Constant wavelength intervals',
'Wavelength $\mu$m')
p.plot(2,(wnc[1:]+wnc[:-1])/2, np.diff(wnc,axis=0), 'Non-constant wavenumber intervals',
'Wavenumber cm$^{-1}$');






## Mixing spectral data from different sources

When working with mixed domains, e.g., Modtran transmittance data in wavenumber domain and sensor spectral response in wavelength domain, you have to resample from the one to the other, to at least get samples at the same spectral intercepts. So, work in either domain, but be sure that the sampled data all relate to the same fundamental spectral baseline. The pyradi library provides functions to resample data when reading in the data from file, with minimal coding overhead (see ryfiles.loadColumnTextFile).

## Performing spectral integrals

A repeating pattern in computational radiometry is the spectral integral: integrating the sum of a spectral variable between an upper and a lower limit. For example, this calculation is required to determine the total inband radiance over a wide spectral range. We use the mathematical integral notation to indicate such integrals, but in practice, these are calculated using numerical summation. $$L = \int_{\lambda_0}^{\lambda_1} L_\lambda d\lambda = \sum_{i=0}^{N}L(i)\Delta\lambda(i)$$ where the sample $i=0$ corresponds with the sample at $\lambda_0$ and the sample $i=N$ corresponds with the sample at $\lambda_1$ The sum above is a general statement and does not make any assumption about the spectral domain sampling interval. If the spectral domain sampling interval has fixed increments $\Delta\lambda(i)=\Delta\lambda$ then $$L = \int_{\lambda_0}^{\lambda_1} L_\lambda d\lambda = \Delta\lambda\sum_{i=0}^{N}L(i)$$

In pyradi the spectral integrals are typically (and naively) calculated using the numpy.trapz trapezoidal integration function. This function has the signature
numpy.trapz(y, x=None, dx=1.0, axis=-1)
where y is the value to be integrated (i.e., spectral density $L_\lambda$), x is the domain where the value y is sampled(i.e., wavelength $\lambda$), dx is the constant $\Delta\lambda$ above, and axis is the axis along which integration is done.

If the spectral domain is sampled at regular intervals the x array is not not required and unless dx is specified, it is assumed that dx=$\Delta\lambda$=1. In this case the spectral sample interval must be accounted for, either by (1) pre-multiplication, (2) by stating the value for dx or by (3) supplying the value for x:



In [13]:

wl = np.linspace(1, 14, 1001)#.reshape(-1,1) # wavelength
temperature = 1000.0 # K
radiancewl = ryplanck.planck(wl,temperature, 'el') / np.pi
print('Integral(1) = {} W/(m^2.sr)'.format( (wl[1]-wl[0]) * np.trapz(radiancewl)) )
print('Integral(2) = {} W/(m^2.sr)'.format( np.trapz(radiancewl, dx=(wl[1]-wl[0]))) )
print('Integral(3) = {} W/(m^2.sr)'.format( np.trapz(radiancewl,wl)) )




Integral(1) = 17373.03021854981 W/(m^2.sr)
Integral(2) = 17373.030218549815 W/(m^2.sr)
Integral(3) = 17373.030218549946 W/(m^2.sr)



If the spectral interval is not constant, i.e, wavelength domain converted to wavenumber domain, the dx or pre-multiplication approach is no longer valid. The only possibility to calculate the integral is to give the value of x, with its non-constant intervals. Continuing from the previous example, convert the wavelength domain to the wavenumber domain, calculate the spectral density in wavenumber units W/(m$^2$.sr.cm$^{-1}$) and calculate the integral using the (non-constant-interval) wavenumber domain as x:



In [14]:

wn = ryutils.convertSpectralDomain(wl, type='ln') # wavenumber
radiancewn = ryplanck.planck(wn, temperature, 'en') / np.pi
print('Integral(wn) = {} W/(m^2.sr)'.format( np.trapz(radiancewn, wn )))




Integral(wn) = -17373.300817295316 W/(m^2.sr)



What happened here? The integral appears to have (practically) the same absolute values, but the sign is wrong!

Recall from Section Spectral-variables above, that $d \tilde{\nu} = -d\lambda \,10^4/\lambda^2=-d\lambda\, \tilde{\nu}^2 /10^4$. We can see this by taking the first few elements of the numpy.diff along the spectral domains. Notice also that the wavelength interval is fixed, but the wavenumber interval varies.



In [15]:

print('First 5 intervals in wl: {}'.format(np.diff(wl)[:5]))
print('First 5 intervals in wn: {}'.format(np.diff(wn)[:5]))
#now calculate the intervals using the equations above:
print('First 5 intervals in wn: {}'.format( - np.diff(wl)[:5] * 1.0e4 / ((wl[:5] + wl[1:6])/2. )**2 ))
print('First 5 intervals in wn: {}'.format( - np.diff(wl)[:5] * ((wn[:5] + wn[1:6])/2. )**2 / 1.0e4))




First 5 intervals in wl: [ 0.013  0.013  0.013  0.013  0.013]
First 5 intervals in wn: [-128.33168806 -125.07961799 -121.94961792 -118.93565398 -116.03206055]
First 5 intervals in wn: [-128.32633585 -125.0745336  -121.94478481 -118.93105681 -116.0276851 ]
First 5 intervals in wn: [-128.33704049 -125.08470258 -121.95445123 -118.94025133 -116.03643617]



The intervals are correct in absolute magnitude (corresponding to the spectral samples), just opposite in sign. To fix this, just take the negative or the absolute value as in



In [16]:

print('Integral(wn) = {} W/(m^2.sr)'.format( - np.trapz(radiancewn, wn )))
print('Integral(wn) = {} W/(m^2.sr)'.format( np.abs( np.trapz(radiancewn, wn ))))




Integral(wn) = 17373.300817295316 W/(m^2.sr)
Integral(wn) = 17373.300817295316 W/(m^2.sr)



The other alternative is to flip the spectral domain around (to obtain positive increments) and then do the integral. But there is catch here: the intervals are not constant, so you must flip (numpy.flipud) both the spectral density and spectral domain variables, to ensure matching elements.



In [17]:

#flip only wavenumber but not the wavenumber spectral density
print('Integral(wn) = {} W/(m^2.sr) - wrong because of a mismatch between x and y ordering'.\
#flip both wavenumber and wavenumber spectral density
print('Integral(wn) = {} W/(m^2.sr)'.format( np.trapz(np.flipud(radiancewn), np.flipud(wn) )))




Integral(wn) = 27199.45508877303 W/(m^2.sr) - wrong because of a mismatch between x and y ordering
Integral(wn) = 17373.300817295316 W/(m^2.sr)



## Converting wideband values between radiant and photon rate units

The relationship between the photon rate radiometric units and radiant units are well defined for a single wavelength. It is best to calculate the spectral photon rate units from first principles and then integrate over the spectral band of interest. However on occasion the information may only be available in wideband radiant units, but the quantity must be expressed in photon rate units. For wideband data the question then arises which wavelength should be used in the conversion? This section determined the appropriate wavelength to be used for the conversion for a given spectral band and source spectrum.

The appropriate conversion wavelength is determined at the wavelength where an accurate spectral calculation gives the same answer as an approximate wideband calculation. The experimentation below indicates that the conversion wavelength is slightly longer than the center of the spectral band - photons at longer wavelengths have lower energy than at shorter wavelengths.

Responsivity is given by $$R_e = \frac{\eta q \lambda}{hc}$$ from which quantum efficiency can be determined $$\eta = \frac{R_e hc}{q \lambda}$$



In [18]:

# to investigate the conversion for broadband data between photon rate units and radiant units
import scipy.constants as const

wl = spec[:,0].reshape(-1)
respW = 2. * spec[:,2].reshape(-1)
respQ = respW * const.h * const.c / (wl * 1e-6 * const.e)

wlNan = np.where(respW==0,np.nan,wl)
wlcent = np.mean([np.nanmax(wlNan),  np.nanmin(wlNan)])
print('wl min = {}'.format(np.nanmin(wlNan)))
print('wl mid = {}'.format(wlcent))
print('wl max = {}'.format(np.nanmax(wlNan)))

p = ryplot.Plotter(1,2,2,figsize=(10,6));
p.plot(1, wl,respW, 'Radiant domain','Wavelength $\mu$m','Responsivity A/W');
p.plot(2, wl,respQ,'Wavelength $\mu$m','Quantum efficiency');

q = ryplot.Plotter(2,1,1,figsize=(8,4))
for temp in [250.,300.,350.]:
Ve = respW * ryplanck.planck(wl,temp,'el') / np.pi
Vq = respQ * ryplanck.planck(wl,temp,'ql') / np.pi
Vei = np.trapz(Ve,-wl,axis=0)
Vqi = np.trapz(Vq,-wl,axis=0)

# convert from W to q/s
qW = const.h * const.c /  (wl * 1e-6)
ratio =  ((Vei/qW) / Vqi ) / (np.max(respW)/np.max(respQ))
wlu = np.interp(1.,ratio[::-1],wl[::-1])
q.plot(1, wl,ratio, 'Ratio (normalised responses) as function of conversion wavelength',
'Wavelength $\mu$m', 'Ratio',
label=['{} K, crossover at {:.3f} $\mu$m'.format(temp,wlu)]);

'Wavelength $\mu$m', 'A/(m$^2$.sr.$\mu$m)',
label=['{} K, {:.3e} A/(m$^2$.sr)'.format(temp,Vei)]);
p.plot(4, wl,Vq, 'Quantum domain ',
'Wavelength $\mu$m', 'q/(s.m$^2$.sr.$\mu$m)',
label=['{} K, {:.3e} q/(s.m$^2$.sr)'.format(temp,Vqi)]);
#     p.plot(4, wl,Ve/const.e);




wl min = 3.502627
wl mid = 4.1963013
wl max = 4.8899756



## Effective value normalisation

The effective value of a spectral variable is given by (Sec 7.2.2)

$${F}_{\rm eff} = \frac {\int_{0}^{\infty}{F} {G} d\lambda} {\int_{0}^{\infty}{G} \lambda},$$

where ${F}$ is the variable in question, and ${G}$ is a weighting function. Note that the effective value of $F$ depends on the shapes of both $F$ and $G$; the effective value of $F$ thus calculated therefore applies only to the specific weighting function $G$ - it is incorrect for any other shapes.

To demonstrate the sensitivity of the effective value normalisation to the spectral variability of the input spectra, consider the effective detector responsivity calculated at different blackbody source temperatures from 500 K to 2000 K. The InSb detector is modelled using the simple detector model in ryutils.responsivity. The source radiance is the weighting function $G$, calculated by using ryplanck.planck. The following graphs show the spectral detector responsivity ($F$ in the above equation) and the effective responsivity as calculated for the range of temperatures. The calculation requires the calculation of a separate quotient and integral for each temperature, yet there are no loops! This magic is achieved by constructing two-dimensional arrays with wavelengths along the row (0'th) axis and temperature along the column (1'st) axis. The first of these arrays is returned by the ryplanck.planck function (because a vector of wavelengths and a vector of temperatures are passed to the function). The second array is obtained by using the numpy.meshgrid array (google to see how this function works). The two arrays are used to calculate the numerator of the quotient. The integrals are implemented using the numpy.trapz function, integrating along the wavelength axis=0 direction. This integration returns a vector of values corresponding to the temperature values along the columns axis=1.

Note the very wide variation in effective responsivity when calculated with different source temperatures. This stems from the fact that the responsivity is weighted differently by the different Planck radiance curves resulting from the different source temperatures. The highest effective responsivity occurs near the wavelength where the peak of the Planck radiance matches the peak of the detector responsivity.



In [19]:

wavelength = np.linspace(0.1, 20, 305).reshape(-1, 1)
detector = ryutils.responsivity(wavelength, 6.1, 1, 15, 1.0)
temperatures = np.linspace(500, 2000, 41)
spectralBaseline = ryplanck.planck(wavelength,temperatures,'el')
_, detVar = np.meshgrid(temperatures, detector)
spectralWeighted = spectralBaseline * detVar
effResp = np.trapz(spectralWeighted, wavelength, axis=0 ) / np.trapz(spectralBaseline, wavelength, axis=0)

p = ryplot.Plotter(2,1,3,figsize=(18,6))
p.plot(1,wavelength,detector,'Spectral responsivity','Wavelength $\mu$m', 'Responsivity A/W',pltaxis=[0.5,7,0,0.8])
p.plot(2,temperatures,effResp,'Effective responsivity','Temperature K', 'Responsivity A/W')




Out[19]:

<matplotlib.axes._subplots.AxesSubplot at 0xd9fc518>



Data are easily read in from files. The format most easy to read/write is ASCII coding in a tabular format (also known as a flat file). ASCII tabular files can use commas (csv format) or whitespace (tabs, spaces) to delimit the columns from each other. Python provides specialist CSV readers or just ASCII format readers. The data processing package Pandas also has very powerful data reading and manipulation facilities, especially for data in complex structural formats. In the pyradi context we are more interested in opening files with numpy and this analysis will focus on reading and writing files for Numpy applications.

Numpy (in versions later than 1.6) has powerful functions to read and write data files to and from Numpy arrays. The file formats include ASCII, binary format and compressed binary formats.

ASCII flat files are convenient data stores because these file can be read/written by a text editor and with a little effort with Excel (import data from txt, save as txt). In this example an array is created, written to file and read back from the file. When writing the text file, the string format is set to truncate the number of decimal values - on purpose, because we want to compare the array in memory with the array read in from file. The numpy.loadtxt and numpy.savetxt functions does the work. These functions take as a parameter the filename or file handle (if the file is opened elsewhere). If the filename has the extension .gz, the file is handled as a gzip compressed file.

When reading and writing files, the with context manager is used. The motivation for using with is that if something goes wrong in the opening, reading or writing, the file is closed in an orderly manner. Also, it is more concise and requires less typing, because the with context manager for the file closes the file automatically when leaving the block. Note that the with context manager provides a file handle: the gz extension is lost and the file is not compressed.

In this example the file is written a two-line header, started with the comment symbol. When the file is read again later, the lines with the comment symbol are ignored. The last few lines read only selected columns - useful if you do not need all the data in the file.



In [20]:

arr = np.random.randn(3,4)
print('original data:\n{}'.format(arr))

filename = 'arrfile.txt'
with open(filename, 'wb' ) as filehandle:
np.savetxt(filehandle, arr, fmt='%12.5e',
with open(filename, 'rb' ) as filehandle:
print('read from text file {} data:\n{}'.format(filename, arr2))

#write/read gzip compressed file - does not work with with
filename = 'arrfile.gz'
np.savetxt(filename, arr, fmt='%12.5e',
print('read from text file {} data:\n{}'.format(filename, arrz))
usecols=[0,3]
print('read columns {} from text file {} data:\n{}'.format(usecols,
filename, arrz))




original data:
[[ -9.91893151e-02   4.41916245e-01  -9.85053886e-01  -2.40719391e-01]
[ -1.33119164e+00   5.93467897e-01   1.30289892e-01  -5.33955740e-01]
[ -6.50621499e-02   1.20936124e-03  -1.78552080e-01   1.97165411e-01]]
read from text file arrfile.txt data:
[[ -9.91893000e-02   4.41916000e-01  -9.85054000e-01  -2.40719000e-01]
[ -1.33119000e+00   5.93468000e-01   1.30290000e-01  -5.33956000e-01]
[ -6.50621000e-02   1.20936000e-03  -1.78552000e-01   1.97165000e-01]]
read from text file arrfile.gz data:
[[ -9.91893000e-02   4.41916000e-01  -9.85054000e-01  -2.40719000e-01]
[ -1.33119000e+00   5.93468000e-01   1.30290000e-01  -5.33956000e-01]
[ -6.50621000e-02   1.20936000e-03  -1.78552000e-01   1.97165000e-01]]
read columns [0, 3] from text file arrfile.gz data:
[[-0.0991893 -0.240719 ]
[-1.33119   -0.533956 ]
[-0.0650621  0.197165 ]]



Data can also be written in binary format or compressed binary format. These files can only be read by the appropriate binary reading files, and not with a text editor. As a rule, these files are smaller than text files. The binary file format does not support headers and is written and read as a single block of data.

In this example the compressed file is not much smaller than the uncompressed file, because the data are random numbers - compression does not help much for this type of data.



In [21]:

arr = np.random.randn(3,4)
print('original data:\n{}'.format(arr))

filename = 'arrfile.npy'
np.save(filename, arr)
print('read from binary file {} data:\n{}'.format(filename, arr2))

filename = 'arrfile.npz'
np.savez(filename, arr=arr)
print('read from text file {} data:\n{}'.format(filename, arrz['arr']))

filename = 'arrfilez.npz'
np.savez_compressed(filename, arr=arr)
print('read from text file {} data:\n{}'.format(filename, arrz['arr']))




original data:
[[-0.21861692  0.17650595  0.13094845  0.55736021]
[-1.18510384 -1.97529995 -0.45643309 -0.0590228 ]
[-1.19088657  0.17834814  2.22421946 -0.38280804]]
read from binary file arrfile.npy data:
[[-0.21861692  0.17650595  0.13094845  0.55736021]
[-1.18510384 -1.97529995 -0.45643309 -0.0590228 ]
[-1.19088657  0.17834814  2.22421946 -0.38280804]]
read from text file arrfile.npz data:
[[-0.21861692  0.17650595  0.13094845  0.55736021]
[-1.18510384 -1.97529995 -0.45643309 -0.0590228 ]
[-1.19088657  0.17834814  2.22421946 -0.38280804]]
read from text file arrfilez.npz data:
[[-0.21861692  0.17650595  0.13094845  0.55736021]
[-1.18510384 -1.97529995 -0.45643309 -0.0590228 ]
[-1.19088657  0.17834814  2.22421946 -0.38280804]]



pyradi.ryplot also provides specialist functionality to read files. The use of these files is demonstrated below in the colour example.

1. The function ryfiles.loadColumnTextFile was written to load spectral files in pyradi. The first column contains the spectral domain (e.g., wavelength or wavenumber), and the remaining columns contain the spectral variables (e.g.. transmittance). The spectral domain can be scaled; from nanometre to micormetre. The spectral variables can be interpolated to a new spectral domain. This is often required because the spectral domain values used in a calculation is not the same as the spectral domain used in the data file. In short using pyradi.ryfiles.loadColumnTextFile can shorten the code significantly when reading in spectral data files.

2. It is common practice to document the data in a file with a one-line header containing the names of the columns in the file. The function ryfiles.loadHeaderTextFile reads the column headers from the first line of a text file. The column names in this header file must be comma separated, because some column headings may contain spaces.

3. A function is provided to take an arbitrary string and turn it into a legal filename by removing characters not allowed in filenames. A default set is provided but you may provide your own. See ryfiles.cleanFilename.

4. Two functions are provided to unzip gzip files ryfiles.unzipGZipfile and untar tar files ryfiles.untarTarfile. Finally, a function is provided to download a file from the internet, given the URL ryfiles.downloadFileUrl.

## Spectral convolution

(Sec 6.6) describes how to calculate the signal that a given sensor with response $S$ would receive from a source $L$, through some medium with transmittance $\tau_{01\lambda}$. Suppose that this sensor has a spectral filter $\tau_f$ with a narrow (but nonzero) spectral width. This filter is very narrow compared with its central wavelength, say $\Delta \lambda = 0.01 \lambda_c$. Such radiometers are used to determine the spectral radiance of sources or to measure the spectral transmittance of the atmosphere. The apparent irradiance measured by such a system can be written [from (Eq 6.13)] $$E_{\lambda_c}= k\int_{0}^{\infty} \epsilon_{0\lambda} L_{0\lambda} \,\tau_{a\lambda} \,\tau_{f\lambda} {S}_\lambda \lambda,$$ where $k$ accounts for the geometrical factors such as source area, orientation, and distance. This equation can be written as $$E_{\lambda_c}= k\int_{\lambda_c - \frac{\Delta\lambda}{2}}^{\lambda_c+\frac{\Delta\lambda}{2}} \epsilon_{0\lambda} L_{0\lambda} \,\tau_{a\lambda} \,\tau_{f\lambda} { S}_\lambda \lambda.$$ By change of variable $\lambda = \lambda_c-x$, $$E_{\lambda_c}= k\int_{-\frac{\Delta\lambda}{2}}^{+\frac{\Delta\lambda}{2}} \epsilon_{0x} L_{0x} \,\tau_{ax} \,\tau_{f(\lambda_c-x)} { S}_x x.$$ These equations show very clearly that the irradiance measured with the filter centred around wavelength $\lambda_c$ includes source energy from $\lambda_c - \frac{\Delta\lambda}{2}$ to $\lambda_c+\frac{\Delta\lambda}{2}$. Apart from the spectral selection, the filter has an additional effect by smoothing the spectrum being observed because the filter has a nonzero spectral width.

The equation above is called a convolution integral because it describes the convolution between the product $(\epsilon_{0\lambda_c}L_{0\lambda_c} \,\tau_{a\lambda_c} \, { S}_{\lambda_c})$ and $\tau_{f}$. In linear systems terminology, the observed spectral source radiance is being convolved with the filter spectral transmittance. To investigate the effects of this convolution consider the two cases: (1) the observed source has little variation over the filter passband, and (2) the observed source varies significantly over the filter passband.

If the product $(\epsilon_{0\lambda_c} L_{0\lambda_c} \,\tau_{a\lambda_c} \, { S}_{\lambda_c})$ is more or less constant over the filter passband $\Delta\lambda$, the equation can be written to show that the convolution has little effect other than some insignificant amount of smoothing: $$# E_{\lambda_c} k\,(\epsilon_{0\lambdac} L{0\lambdac} \,\tau{a\lambdac} \, { S}{\lambdac}) \int{-\frac{\Delta\lambda}{2}}^{+\frac{\Delta\lambda}{2}} \tau_{f(\lambdac-x)} x \approx \, k\,\epsilon{0\lambdac} L{0\lambdac} \,\tau{a\lambdac} \, { S}{\lambdac}\, \tau{f}\,\Delta x.$$ If the product $(\epsilon\lambda L{\lambda}\;\tau{a\lambda}\;{ S}\lambda)$ varies significantly over the filter passband $\Delta\lambda$, the convolution equation cannot be simplified. In this case, the convolution attenuates and smears out the finer detail in the spectral information. If the actual spectral line is very narrow the measured line will approximate the filter resolution and will be totally erroneous unless the filter spectral smear effect is compensated by deconvolution.

This theory may be intellectually satisfying, but it is much nicer to see the convolution principle in graphical form with real-world information. First investigate the principle of convolution by creating a synthetic signal by placing a number of impulses at varying distances from each other. Then convolve this signal by a wider top-hat (square) function. Note how the top-hat function 'smears out' the signal. The pyradi.ryutils.convolve function is used to do the convolution.



In [22]:

samplingresolution=0.5
wavenum=np.linspace(0, 100, int(100/samplingresolution))
inspectral=np.zeros(wavenum.shape)
inspectral[int(10/samplingresolution)] = 1
inspectral[int(11/samplingresolution)] = 1
inspectral[int(30/samplingresolution)] = 1
inspectral[int(45/samplingresolution)] = 1
inspectral[int(55/samplingresolution)] = 1
inspectral[int(70/samplingresolution)] = 1
inspectral[int(75/samplingresolution)] = 1
inwinwidth=1
outwinwidth=5
outspectral,  windowfn = ryutils.convolve(inspectral, samplingresolution,  inwinwidth,  outwinwidth)
convplot = ryplot.Plotter(1, 1, 1, figsize=(8,4))
convplot.plot(1, wavenum, inspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['r-'],label=['Input'],legendAlpha=0.5)
convplot.plot(1, wavenum, outspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['g-'],label=['Output'],legendAlpha=0.5);






### Savitzky-Golay filter

The Savitzky-Golay filter smooths signals by by fitting successive sub-sets of adjacent data points with a low-degree polynomial by the method of linear least squares.

A Python implementation of this filter is given here. A one-dimensional version of the filter is also included in pyradi.ryutils.SavitzkyGolay.

The example below shows a considerable ringing on the spiky input signal. Clearly this filter is not appropriate for this particular input signal. Another example (atmospheric transmittance) is shown below which better demonstrates the Savitzky-Golay filter performance.



In [23]:

outspectral = ryutils.savitzkyGolay1D(inspectral, window_size=21, order=3, deriv=0, rate=1)

convplot = ryplot.Plotter(1, 1, 1, figsize=(8,4))
convplot.plot(1, wavenum, inspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['r-'],label=['Input'],legendAlpha=0.5)
convplot.plot(1, wavenum, outspectral, "Convolution Test", r'Wavenumber cm$^{-1}$',\
r'Signal', ['g-'],label=['Output'],legendAlpha=0.5);






The next example uses spectral Bunsen flame measurement data. The data is read in from the pyradi Google Code site, so you must be online to run this cell. The data is read using the urllib2 library, given the URL, the urllib2 returns a file handle that is used by Numpy to read the data. The spectral width for the Bunsen flame measurement was 4 cm$^{-1}$ and the signal was sampled at 2 cm$^{-1}$. The Bunsen data file format is: first column contains the wavenumber and the second column contains the raw signal (proportional to radiance). Modtran transmittance data were calculated for a 5 m path length and a spectral width of 1 cm$^{-1}$ and sampling of 1 cm$^{-1}$. The transmittance file format is: first column contains the wavenumber and the second column contains the transmittance. The atmospheric transmittance is convolved to resolutions of 4 cm$^{-1}$ and 40 cm$^{-1}$. The spectral baselines for the Bunsen data and Modtran data are not the same, so the Bunsen data is re-sampled at the Modtran spectral values. All the data is then plotted for visual inspection. The convolution smoothing effect is clearly visible in the graphs.



In [24]:

from scipy.interpolate import  interp1d

wavenum =  tauAtmo[:, 0]
tauA = tauAtmo[:, 1]

# convolve transmittance from 1cm-1 to 4 cm-1
tauAtmo4,  windowfn = ryutils.convolve(tauA, 1, 1, 4)
# convolve transmittance from 1cm-1 to 40 cm-1
tauAtmo40,  windowfn = ryutils.convolve(tauA, 1, 1, 40)

#interpolate bunsen spectrum to atmo sampling
#first construct the interpolating function, using bunsen
#then call the function on atmo intervals
bunsen = bunInterp1(wavenum)

atmoplot = tauA.copy()
atmoplot =  np.vstack((atmoplot, tauAtmo4, tauAtmo40))
convplot02 = ryplot.Plotter(1, 1, 2,figsize=(20,5))
convplot02.plot(1, wavenum, atmoplot.T, "Atmospheric Transmittance", r'Wavenumber cm$^{-1}$',
r'Transmittance', ['r-', 'g-','b-'],label=['1 cm-1', '4 cm-1', '40 cm-1' ],legendAlpha=0.5);
convplot02.plot(2, wavenum, atmoplot.T, "Atmospheric Transmittance", r'Wavenumber cm$^{-1}$',
r'Transmittance', ['r-', 'g-','b-'],label=['1 cm-1', '4 cm-1', '40 cm-1' ],legendAlpha=0.5,
pltaxis=[3700, 3900, 0, 1]);






The next step is to attempt to correct for the atmospheric transmittance present during the measurement. In this particular case the simple approach of dividing the radiance signal by the transmittance is attempted. The case on the left side is for the Bunsen data at 4 cm$^{-1}$ resolution and the atmospheric transmittance (mismatched at) at 1 cm$^{-1}$. The case on the right side is for the Bunsen data at 4 cm$^{-1}$ resolution and the atmospheric transmittance (well matched at) at 4 cm$^{-1}$. The results from the mismatched calculation is clearly much more noisy than the results for the well-matched calculation. The results also show that it is impossible to reconstruct the signal in spectral regions where there is zero signal - a more sophisticated model-based method must be followed, such as explained in (Sec 8.4).



In [25]:

bunsenPlt = ryplot.Plotter(1,3, 2, figsize=(15,7))
bunsenPlt.plot(1, wavenum, bunsen, "Bunsen Flame Measurement 4 cm-1", r'',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(2, wavenum, bunsen, "Bunsen Flame Measurement 4 cm-1", r'',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(3, wavenum, tauA, "Atmospheric Transmittance 1 cm-1", r'',
r'Transmittance', ['r-'])
bunsenPlt.plot(4, wavenum, tauAtmo4, "Atmospheric Transmittance 4 cm-1", r'',
r'Transmittance', ['r-'])
bunsenPlt.plot(5, wavenum, bunsen/tauA, "Atmospheric-corrected Bunsen Flame Measurement 1 cm-1", r'Wavenumber cm$^{-1}$',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5])
bunsenPlt.plot(6, wavenum, bunsen/tauAtmo4, "Atmospheric-corrected Bunsen Flame Measurement 4 cm-1", r'Wavenumber cm$^{-1}$',
r'Signal', ['r-'], pltaxis =[2000, 4000, 0,1.5]);