4 Introduction to computational radiometry with pyradi

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.

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 in pyradi

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.

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:

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.

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.

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.

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.

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.

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.

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).

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:

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:

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.

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

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.

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 \begin{equation} R_e = \frac{\eta q \lambda}{hc} \end{equation} from which quantum efficiency can be determined \begin{equation} \eta = \frac{R_e hc}{q \lambda} \end{equation}

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.

File loading and saving

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.

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.

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.

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.

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.

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).

The atmospheric transmittance model is next Savitzky-Golay filtered. This filter appears to be popular in the spectroscopy community, so it should work well for the atmospheric transmittance data used here. It indeed does provide nicely filtered spectra, as seen below. More work is required to accurately translate the SG filter parameters to spectral widths.

Colour calculations

This is a very high-level introduction to colour theory. For more detail see
https://github.com/colour-science/colour
http://colour-science.org/tutorial.php
http://nbviewer.ipython.org/github/colour-science/colour-ipython/blob/master/notebooks/colour.ipynb
http://nbviewer.ipython.org/github/colour-science/colour-ipython/blob/master/notebooks/colorimetry/spectrum.ipynb

(Sec 2.10.5) Colour is an elusive property - different people perceive colour differently, and the apparent colour of an object depends on the illuminance spectrum. The colour tristimulus have been standardised by the CIE. It is interesting to note that the red colour stimilus has a significant component towards the blue part of the spectrum. The description in this section is very brief, for more details see (Sec 2.10). The example given here investigates the human perception of colour, relative to the CIE 1931 standard.

This section explores the subtleties of colour as an application of normalising and radiometric concepts and colour coordinates. Four sources are considered: the first light source is a 'daylight' fluorescent light source, the second source is the sun modelled as a thermal radiator at 5900~K, the third source is an incandescent light globe at a temperature of 2850~K, and the fourth source is a low-pressure sodium lamp.

The samples illuminated by these sources are a red tomato, lettuce, a yellow prune, a dark-green leaf, a blue Nitrile (latex-like) surgical glove, and standard white printing paper. The samples' diffuse reflection spectra were measured with a spectroradiometer, illuminating the sample with a bright light at short distance. The fruit samples all demonstrated considerable light propagation deeper into the fruit. The blue glove was located on a Spectralon white reference (note the considerable 'white' reflectance beyond 0.55~$\mu$m).

Use a local copy of the data file, if the tgz file is available in the same directory where the notebook is located. Otherwise download the 28 KB file from the internet. Once availabe, unzip and untar the data files to support rest of the calculation.

The data is read with ryfiles.loadColumnTextFile. This function loads column data from a text file, scaling and interpolating the read-in data, according to user specification. The first 0'th column has special significance: it is considered the abscissa (x-values) of the data set, while the remaining columns are any number of ordinate (y-value) vectors. The user passes a list of columns to be read (default is [1]) - only these columns are read, processed and returned when the function exits. The user also passes an abscissa vector to which the input data is interpolated and then subsequently amplitude scaled or normalised.

Next the source spectra are loaded and/or calculated. The following code builds a 2-D array with the source radiance values, where each column represents a different source. Wavelength extends along rows. Spectral interval for all source spectra is the same, which is the wavelength array calculated in the previous cell. Blackbody radiance spectra are calculated at the required wavelength intervals for the required temperatures. All data read from files are interpolated to the required wavelength intervals

The spectra from the four sources are stacked horizontally by using numpy.hstack. The resulting array varies with wavelength along the rows (axis=0), with the four source spectra are given in the four columns.

Next load the sample reflectance data. The data file for this data has the wavelength ($\mu$m) values along the first column, with many columns following to represent the different samples. Only a few of the samples (samplesSelect = [1,2,3,8,10,11]) are loaded here.

After reading in the numeric values, the headers in the first row is read in. These headers are used in labels in the graphs.

Now calculate the CIE x,y coordinates for each combination of source and sample. The x,y coordinates indicate the colour of the sample under the specific source illumination. This cell uses the rychroma.chromaticityforSpectralL function to do the coordinate calculation from the spectra. For equations and more details see (Sec 2.10).

We have a number of samples and a number of sources, essentially these two are orthogonal to each other. The following calculation is therefore done in two loops and the results written to two-dimensional arrays for storage and subsequent plotting. The two arrays are pre-allocated and filled with zeros before the calculation.

This cell calculates the horseshoe-shape locus of monochromatic light. The concept is really quite simple: create a monochromatic spectrum (only one wavelength) and then use rychroma.chromaticityforSpectralL to calculate the x,y coordinates for that wavelength.

Finally plot the colour coordinates for each sample-source combination on the diagram, together with the monochromatic horseshoe. the graph is plotted trice; once to get the big picture and then two zoomed in on the detail.

It is quite evident that for the two 'daylight' sources (the sun 5900 K and the fluorescent source) the sample colour coordinates are widely scattered over the coordinate space. When the same samples are viewed under 2850 K source illumination (incandescent lamp) the colours tend to cluster towards the yellow - this matches the human experience that incandescent lamps produce yellow colour photographs. Finally, when viewed under sodium illumination the sample colour coordinates converge on the colour of the lamp, meaning that the samples all appear yellow in colour (589 nm) and only differ in shade.

The CIE 1931 chromaticity diagram in colour is shown below. It was obtained from https://en.wikipedia.org/wiki/CIE_1931_color_space

Solid Angles and Spatial Integrals

Computational radiometry has a strong element of geometry: sources and receivers have size, shape and orientation and are related to each over in some geometrical relationship. The calculation of solid angle is an essential task in most radiometric analyses. The calculation of a solid angle requires integration across a spatial extent as described in (Sec 2.5.4). Array processing in Python condenses and simplifies most two-dimensional calculations significantly.

The solid angle of a rectangular flat surface, with dimensions $W$ and $D$, as seen from a reference point centered above the surface, is determined by the integral of the projected area of a small elemental area $\cos\theta\,dd\,dw$ across the full size of the surface: $$ \omega_{\rm s}=\int_W\int_D\frac{dw\,dd\,\cos^{n-2}\theta}{R^2} $$ $$ \omega_{\rm s}=\int_W\int_D\frac{dw\,dd\,\cos^n\theta}{H^2} $$ $$ \omega_{\rm s}=\int_W\int_D\frac{dw\,dd\,}{H^2}\left(\frac{H}{R}\right)^n $$ $$\omega_{\rm s}=\int_W\int_D\frac{dw\,dd\,}{H^2}\left(\frac{H}{\sqrt{w^2+d^2+H^2}}\right)^n, $$ where $H$ is the reference point height above the surface, and $n=3$ for the geometrical solid angle and $n=4$ for the projected solid angle. The integral is performed along the $W$ and $D$ dimensions with increments of $dw$ and $dd$. The slant range between the reference point and the elemental area $dd\times dw$ is $R=H/\cos\theta$.

This example calculates the geometric and projected solid angles of a 5 m $\times$ 5 m floor from a sensor centered above the floor at a height of 3.5 m. The solid angles thus calculated are compared to the (incorrect) $W\times D/R^2$ approximation of (5 m $\times$ 5 m)/(3.5 m $\times$ 3.5 m). We investigate the accuracy of the answer for a variable number of samples. It seems you need at least 500**2 samples.

Experiment a little with the meshgrid and its results. What do these variables look like graphically? Calculate and plot, for each tile: (1) the distance from the reference point, (2) the solid angle of the angle to the reference point, (3) the geometric solid angle of each tile, and (4) the projected solid angle of each tile. Along the way plot the values of x and y as well.

Compare the projected and geometric solid angles. Note that directly under the reference point the geometric and project solid angles are equal (because $\cos(0)=1$), but at the corners the cosine weighting function starts to weigh the the projected solid angle down.

Another approach, provided by a student, to do the integral as two separate summations. The first summation is in sigmaGSD = sigmaGSD + .... and the second summation is in GeoSA = sum(sigmaGSD). This approach is perhaps more intuitive to some, but requires a for loop [which is perfectly legal, but superfluous in my personal programming mindset].

Python and module versions, and dates