Introduction

In this vignette we will walk through the different steps used in the conversion of raw counts into physical quantities. Different protocols can be used, and not all those implemented are exemplified. The facilities in the package also allow users to implement their own measurement and correction protocols. It must be, however, always kept in mind that protocols and calibrations go hand in hand: For each different measuring protocol or new raw conversion algorithm new calibration multipliers will be required. In the case of new measurement protocols new calibration data and validation measurements will be needed. However, for new raw data conversion protocols recalculation of new calibration multipliers from existing raw-counts and calibration-lamp data may be enough. In this last case existing validation measurements against other instruments can be reused if available.


Code examples

All examples in vignettes use data objects included in the package.

Irradiance

We will step through the different stages of the algorithm used, plotting the spectra at each step along the way. We use raw counts data, included in the package, to make building this vignette easier. At the same time we show how raw data acquired outside R can be processed with the package.

Spectral irradiance from RAW detector counts

We start by loading the R packages we will use.

library(ggplot2)
library(ggspectra)
## Loading required package: photobiology
## Loading required package: tibble
## For news and documentation about the R for Photobiology packages, please, visit http://www.r4photobiology.info/
library(photobiology)
library(photobiologyWavebands)
## For news about 'photobiologyWavebands', please, see http://www.r4photobiology.info/
## For on-line documentation see http://docs.r4photobiology.info/photobiologyWavebands/
library(ooacquire)
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
library(magrittr)

Household white LED bulb

We read RAW-detector-counts data from the measurement of irradiance under a white LED light source—Osram LED Star classic A 60 E27. Spectral data were acquired with an Ocean Optics Maya 2000 Pro spectrometer, using function acquire_irrad_interactive() from this package. This function implements several measurement protocols from which the user can chose. However, RAW counts data returned as a collection of RAW spectra, contains also metadata including a descriptor of the instrument and a descriptor of the instrument settings used.

The data contains a measurement of the light source and a reference dark measurement.

names(white_LED_2min.raw_mspct)
## [1] "light" "dark"

The spectrum for the light measurement contains multiple columns.

names(white_LED_2min.raw_mspct[["light"]])
## [1] "w.length" "counts_1" "counts_2" "counts_3"

Summary in addition to displaying the summary for the columns, displays the most important metadata attributes.

summary(white_LED_2min.raw_mspct[["light"]])
## Summary of object: raw_spct [2,068 x 4]
## Wavelength range 187.82 to 1117.14 nm, step 0.41 to 0.48 nm
## Label: LED_lamp04_long 
## Measured on 2017-04-12 16:34:00 UTC 
## Data acquired with 'MayaPro2000' s.n. MAYP11278
## grating 'HC1', slit '010s'
## integ. time (s): 1.52, 6, 7.06
## total time (s): 120, 120, 120
## counts @ peak (% of max): 94.2
##     w.length         counts_1        counts_2        counts_3    
##  Min.   : 187.8   Min.   : 2148   Min.   : 2026   Min.   : 1997  
##  1st Qu.: 431.7   1st Qu.: 4453   1st Qu.:11141   1st Qu.:12839  
##  Median : 668.7   Median : 4663   Median :11968   Median :13814  
##  Mean   : 663.3   Mean   :12652   Mean   :26127   Mean   :28165  
##  3rd Qu.: 897.6   3rd Qu.:13436   3rd Qu.:45138   3rd Qu.:52280  
##  Max.   :1117.1   Max.   :59705   Max.   :64000   Max.   :64000

Here we extract a single metadata item from one of the spectra.

attr(white_LED_2min.raw_mspct[["light"]], which = "instr.desc")$spectrometer.name
## [1] "MayaPro2000"

To see every single data and metadata item stored we can use str() as shown below, but not run given the length of the output.

# not run
summary(white_LED_2min.raw_mspct[["light"]])

In normal use, we calculate irradiance from a set of raw-counts spectral data using the high-level function s_irrad_corrected(). We pass as first argument the object containing the RAW counts and corresponding metadata and the correction method to be used in the conversion of the RAW counts into irradiance.

The R object returned by s_irrad_corrected() belongs to class "source_spct" (defined in package ‘photobiology’) for which summary and plotting methods are available. We first have a quick look at this object. In the output bellow we will notice that the six RAW measurements, three in the light and three in the dark have been used to calculate a single estimate of spectral irradiance. We can also see that the range of wavelengths is narrower as only wavelengths for which calibration data exists have been retained.

summary(LED_lamp_recalc.spct)
## Summary of object: source_spct [1,421 x 2]
## Wavelength range 251.16 to 898.81 nm, step 0.43 to 0.48 nm
## Label: LED_lamp04_long 
## Measured on 2017-04-12 16:34:00 UTC 
## Time unit 1s
## 
##     w.length       s.e.irrad         
##  Min.   :251.2   Min.   :-2.878e-03  
##  1st Qu.:418.2   1st Qu.:-5.892e-05  
##  Median :582.2   Median : 3.648e-03  
##  Mean   :579.8   Mean   : 2.116e-02  
##  3rd Qu.:742.5   3rd Qu.: 3.406e-02  
##  Max.   :898.8   Max.   : 9.657e-02

By plotting the spectral irradiance data, either in whole or a range of wavelengths we can better see the shape of the emission spectrum of the lamp, and in this example also obtain summaries for the irradiance in different wavebands.

plot(LED_lamp_recalc.spct)

What we show bellow are the individual computation steps that the RAW spectral data have gone through to obtain the spectral irradiance spectrum shown in the figure above. We next print the collection of spectra containing raw counts data. The protocol used consisted in one light measurement and one dark measurement, each with bracketing of the integration time.

We have described above the structure of the object containing the RAW counts data. As we saw above, measurements consisted in three measurements using three different integration times. In addition to different integration times, the number of “scans” was adjusted so that the three spectra consist in the average of scans covering the same total length of time.

Plotting the three bracketed spectra shows the effect of using the optimal integration time plus two longer integration times. To understand the figure one needs to be aware that the spectrometer detector saturates at approximately 64000 counts. In other words if more photons imping the detector cell/pixel than needed for 64000 counts, the reading remains at 64000 counts. This is usually called “signal clipping”.

plot(white_LED_2min.raw_mspct[["light"]])

In the figure above, we can see that in addition to clipping, increasing the integration time, at least in some spectrometers significantly increases the “noise floor” or dark readings. This is even clearly seen in the measurements done in darkness. In this spectra, we can easily see some hot pixels in the detector (bad pixels).

plot(white_LED_2min.raw_mspct[["dark"]])

The first step in the processing is removal of bad pixels, as recorded in the calibration data. What this step achieves can be clearly seen by comparing the plots before and after this step.

for (m in names(white_LED_2min.raw_mspct)) {
  white_LED_2min.raw_mspct[[m]] <-
    skip_bad_pixs(white_LED_2min.raw_mspct[[m]])
  print(plot(white_LED_2min.raw_mspct[[m]]) + ggtitle(m))
}

The second step is to replace saturated (clipped) pixel data, with the missing data marked NA. As NA values are not plotted pixels exactly equal to the maximum possible reading disappear from the plot.

for (m in names(white_LED_2min.raw_mspct)) {
  white_LED_2min.raw_mspct[[m]] <-
    trim_counts(white_LED_2min.raw_mspct[[m]])
  print(plot(white_LED_2min.raw_mspct[[m]]) + ggtitle(m))
}
## Warning: Removed 844 rows containing non-finite values (stat_peaks).

An important consideration is that when a pixel well fills some of the excessive charge migrates to nearby pixels. How many nearby pixels are affected depends on the detector type, but it is at most a few tens of pixels. So in this third step, pixels neighbouring those set to NA in the second step are also set to NA.

for (m in names(white_LED_2min.raw_mspct)) {
  white_LED_2min.raw_mspct[[m]] <-
    bleed_nas(white_LED_2min.raw_mspct[[m]])
  print(plot(white_LED_2min.raw_mspct[[m]]) + ggtitle(m))
}
## Warning: Removed 904 rows containing non-finite values (stat_peaks).

The response of array detectors is not linear in response to the number of photons received. Most precisely when the number of counts gets near the maximum value, the sensitivity to additional photons slightly decreases. This is a results of how “full” the sensor wells are, irrespective of the source of the charge. Consequently this correction should be applied before subtraction of the dark reading. The fourth step is to apply a linearisation function, in most cases supplied with the instrument and possibly stored in the instrument firmware. This function is part of the calibration data for the instrument that is stored as part of the metadata during acquisition of the spectra.

for (m in names(white_LED_2min.raw_mspct)) {
  white_LED_2min.raw_mspct[[m]] <-
    linearize_counts(white_LED_2min.raw_mspct[[m]])
  print(plot(white_LED_2min.raw_mspct[[m]], ylim = c(NA, 65000)) + ggtitle(m))
}
## Warning: Removed 904 rows containing non-finite values (stat_peaks).

The first few pixels in the detector of most array spectrometers are covered and never exposed to light photons. These can be used as a dark reference to correct for thermal noise. This step may appear redundant, but might help in cases when instrument temperature is not the same during the light and dark measurements. An alternative is to use pixels known to be not exposed to radiation from the light source, but within the useful measurement range of the instrument. A white LED lamp is not expected to emit UV-C, so wavelengths between 250 nm and 290 nm can used. Consequently, as the fifth step we remove an estimate of dark signal based on these “non-excited pixels”. As seen in the plots, this step effectively removes the overall dark signal from both light and dark measurements. We can, however, also see that there is residual dark noise remaining at pixel level, not all pixels have the same noise floor. Part of this variation is systematic, and some may be random. As we used a total measurement time of 120 s much of the random noise must have cancelled out.

for (m in names(white_LED_2min.raw_mspct)) {
  white_LED_2min.raw_mspct[[m]] <-
    fshift(white_LED_2min.raw_mspct[[m]], range = c(250,290))
  print(plot(white_LED_2min.raw_mspct[[m]], ylim = c(NA, 65000)) + ggtitle(m))
}
## Warning: Removed 904 rows containing non-finite values (stat_peaks).

At this point we have “clean” RAW counts data. The sixth step is to convert these raw counts into counts per second. As the integration time for each spectrum is stored together with the RAW count data, the function call is simple. After this step, the data acquired using different integration times are expressed in the same units. We take advantage of this to plot using different colours the data from the three different “bracketed” integration times. Blue corresponds to the shortest and optimal integration time, and green to the longest one used. Here we over-plot the three separate spectra.

## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1565.40..8775.73]

## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1565.40..8775.73]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1586.03..6826.44]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1611.57..1000.00]

## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1611.57..1000.00]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1613.97..1000.00]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1611.57..1000.00]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1613.97..1000.00]
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1630.93..1000.00]

The seventh step splicing of the light and dark three counts-per-second spectra each into a single combined spectrum.

for (m in names(white_LED_2min.cps_mspct)) {
  white_LED_2min.cps_mspct[[m]] <-
    merge_cps(white_LED_2min.cps_mspct[[m]])
  print(plot(white_LED_2min.cps_mspct[[m]], ylim = c(NA, 40000)) + ggtitle(m))
}
## Warning in range_check(x, cps.cols): Possible off-range cps values
## [-1630.93..1000.00]

In the eigth step we subtract the dark signal from the light measurement. This reduces the remaining noise.

At this point we have a clean counts-per-second spectrum. The ninth step is to apply the calibration multipliers and any additional instrument specific correction to obtain spectral irradiance in a "source_spct" object.

In an optional, final tenth step we apply smoothing to remove some of the noise.

## 637 possibly 'bad' values in smoothed spectral irradiance
plot(white_LED_2min.spct)

With an LED source, at least in the case of our Maya 2000 Pro, we do not need to worry about stray light. However, when we measure light sources with a significant emission in the infra-red above 1000 nm, stray light becomes a problem. In addition sources with narrow peaks of emission may benefit from deconvolution of a measured slit function.

Tungsten halogen lamp

For our instrument, stray light originates mostly in the infra-red and it can be a serious problem when trying to precisely measuring sunlight or incandescent lamps. For this example we use data acquired using a different protocol, which includes three measurements, being the additional one, a measurement through a filter used to estimate stray light in the UV-region. In this case a UV-cut and IR-pass filter (a piece of clear polycarbonate).

The data contains a measurement of the light source directly and through a filter and a reference dark measurement.

names(halogen.raw_mspct)
## [1] "light"  "filter" "dark"

The spectrum for the light measurement contains two columns with RAW-counts data.

names(halogen.raw_mspct[["light"]])
## [1] "w.length" "counts_1" "counts_2"

Summary in addition to displaying the summary for the columns, displays the most important metadata attributes.

summary(halogen.raw_mspct[["light"]])
## Summary of object: raw_spct [2,068 x 3]
## Wavelength range 187.82 to 1117.14 nm, step 0.41 to 0.48 nm
## Label: Halogen 
## Measured on 2017-03-28 11:24:35 UTC 
## Data acquired with 'MayaPro2000' s.n. MAYP11278
## grating 'HC1', slit '010s'
## integ. time (s): 1.86, 7.2
## total time (s): 5.58, 7.2
## counts @ peak (% of max): 94.4
##     w.length         counts_1        counts_2    
##  Min.   : 187.8   Min.   : 2121   Min.   : 1950  
##  1st Qu.: 431.7   1st Qu.: 5316   1st Qu.:14164  
##  Median : 668.7   Median :22588   Median :64000  
##  Mean   : 663.3   Mean   :25981   Mean   :44225  
##  3rd Qu.: 897.6   3rd Qu.:44996   3rd Qu.:64000  
##  Max.   :1117.1   Max.   :60728   Max.   :64000

As above for the LED lamp, we first calculate spectral irradiance from a set of raw-counts spectral data using the high-level function s_irrad_corrected().

halogen.spct <-
  s_irrad_corrected(halogen.raw_mspct, correction.method= MAYP11278_ylianttila.mthd)
plot(halogen.spct)

In the example above we called the functions used for each of the steps in the computation individually, saving the intermediate results so as to be able to show the partly processed data at each step. The functions, however, support the use of pipes as they all have as their first parameter the one accepting the R object returned by the previous stage.

We first use a “pipe” to apply the same initial processing steps as for the LED bulb data in the previous section. The main difference is that we use as internal instrument dark reference those pixels that never are exposed to radiation. For our instrument they correspond to wavelengths 187.82 nm to 189.26 nm. We obtain this time three spectra containing counts-per-second data.

halogen.cps_mspct <- cps_mspct()
for (m in names(halogen.raw_mspct)) {
  halogen.raw_mspct[[m]] %>%
    skip_bad_pixs() %>%
    trim_counts() %>%
    bleed_nas() %>%
    linearize_counts() %>%
    fshift(range = c(187.82,189.26)) %>%
    raw2cps() %>% 
    merge_cps() -> halogen.cps_mspct[[m]]
}
names(halogen.cps_mspct)
## [1] "light"  "filter" "dark"

We plot the returned spectra, both in full, and the UV region by itself. By careful comparison of these later plots one can see that the signal for filter is larger than for dark. As we know from specifications and measurements that the filter used blocks radiation in this region, the difference is due to stray light (radiation of longer wavelengths being detected as ultraviolet).

for (m in names(halogen.cps_mspct)) {
  print(plot(halogen.cps_mspct[[m]]) + ggtitle(m))
  print(plot(halogen.cps_mspct[[m]], range = c(250, 410)) + ggtitle(m))
}

In the next step we subtract the dark reading from both the light and filter readings, and copy attributes, including instrument settings and calibration data.

## [1] "light"  "filter"
for (m in names(halogen01.cps_mspct)) {
  print(plot(halogen01.cps_mspct[[m]]) + ggtitle(m))
  print(plot(halogen01.cps_mspct[[m]], range = c(250, 410)) + ggtitle(m))
}

As we can see in the plots above, a small amount of stray light is present in both spectra in the UV region. We apply a filter correction using a simple method based on a fixed transmittance value. We set flt.Tfr = 0.9 as these is a good estimate of the transmittance of polycarbonate to the stray light.

## [1] "w.length" "cps"
## [1] "second"
plot(halogen_corrected.cps_spct)

plot(halogen_corrected.cps_spct, range = c(250, 410))

mean(clip_wl(halogen_corrected.cps_spct, range = c(250, 300))[["cps"]])
## [1] 2.747068

The average counts-per-second remaining after correcting for stray light is very small. In the next step we apply the calibration multipliers to obtain spectral irradiance.

## [1] "w.length"  "s.e.irrad"

The spectrum displays some noise at the shortest wavelengths and some interference patterns at the long end. The interference patterns come from the light bulb, but the random noise of increasing amplitude with decreasing wavelength in the UV region is due to the spectrometer.

We can smooth this random noise if desired.

plot(halogen.source_spct)

plot(halogen.source_spct, range = c(250, 410))

We can compute some photon ratios, expressed as \(mmol\, mol^{-1}\) to diagnose whether the stray light is well controlled.

q_ratio(halogen.source_spct, list(UVC(), UVB(), UVA()), PAR()) * 1e3
##  UVC.ISO.tr.lo: PAR(q:q)        UVB.ISO: PAR(q:q)        UVA.ISO: PAR(q:q) 
##                0.1925562                0.3494061                3.2901570 
## attr(,"radiation.unit")
## [1] "q:q ratio"

Smoothing can sometimes help, but it can also introduce bias. It should be used with care and always checking the output. The arguments and method chosen remove random noise and a small “bump” at 320 nm, setting strength to 1 instead of 3 smooths the random noise but not this small peak (not shown). In the infra-red most of the wavy pattern is also removed. So, smoothing can be useful, but it can also remove real features, and one needs to decide if these features are of interest or not, based on other sources of information.

plot(halogen_sm.source_spct)

plot(halogen_sm.source_spct, range = c(250, 410))

How much difference did smoothing do?

q_ratio(halogen.source_spct, list(UVC(), UVB(), UVA()), PAR()) * 1e3
##  UVC.ISO.tr.lo: PAR(q:q)        UVB.ISO: PAR(q:q)        UVA.ISO: PAR(q:q) 
##                0.1925562                0.3494061                3.2901570 
## attr(,"radiation.unit")
## [1] "q:q ratio"

Slit function correction

When measuring spectra containing narrow peaks or steep slopes the shape of the slit function will affect the apparent width of peaks or the apparent steepness of the slope. The slit through which light enters the instrument has a finite width and angle of acceptance for radiation. Consequently, even when measuring true monochromatic radiation from a laser photons will imping on multiple array pixels/wells. The pixel corresponding to the wavelength of the radiation receives the most photons, but the neighbouring pixels receive a decreasing number of photons as the distance from the “correct” target pixel increases. The wider the slit used, the broader the “false” peak observed for monochromatic light. The shape and width of this slit function depends not only on the width of the slit but also on other features of the optical bench of the instrument. In array spectrometers the slit function also depends on wavelength as the length of the path from the grating to the detector is not constant.

However, if the slit function is known, it can be used to remove its influence from a measured spectrum, in simpler words it can be used to partly reconstruct the structure of original light source spectrum. In the case of the two examples above, applying this correction would make little difference. In the case of the solar spectrum and discharge lamps this further step improves the estimates for spectral irradiance.

Too see the effect of this correction we need to look at individual peaks in a spectrum. We use data for a mercury lamp.

The data contains a measurement of the light source and a reference dark measurement.

names(xenon_flash.raw_mspct)
## [1] "light" "dark"

The spectrum for the light measurement contains one column with RAW-counts data as no bracketing was used.

names(xenon_flash.raw_mspct[["light"]])
## [1] "w.length" "counts"

Summary in addition to displaying the summary for the columns, displays the most important metadata attributes.

summary(xenon_flash.raw_mspct[["light"]])
## Summary of object: raw_spct [2,068 x 2]
## Wavelength range 187.82 to 1117.14 nm, step 0.41 to 0.48 nm
## Label: Bare bulb xenon flash 
## Measured on 2018-09-13 14:34:06 UTC 
## Data acquired with 'MayaPro2000' s.n. MAYP11278
## grating 'HC1', slit '010s'
## integ. time (s): 5.08
## total time (s): 5.08
## counts @ peak (% of max): 91.9
##     w.length          counts     
##  Min.   : 187.8   Min.   : 2096  
##  1st Qu.: 431.7   1st Qu.: 8150  
##  Median : 668.7   Median :13992  
##  Mean   : 663.3   Mean   :13795  
##  3rd Qu.: 897.6   3rd Qu.:17637  
##  Max.   :1117.1   Max.   :57992

In the case of these data, the concept of counts-per-second does not apply as the flash discharge is shorter than the integration time, and unknown. The relevant reference is one exposure event and the quantity to estimate is spectral fluence in \(J m^{-2}\).

## NULL
## [1] 1

As above for the LED lamp, we first calculate spectral fluence from a set of raw-counts spectral data using the high-level function s_irrad_corrected().

xenon_flash.spct <-
  s_irrad_corrected(xenon_flash.raw_mspct, correction.method = MAYP11278_ylianttila.mthd)
getTimeUnit(xenon_flash.spct)
## [1] "exposure"
plot(xenon_flash.spct, range = c(315, NA))

xenon_flash.cps_spct <-
  s_irrad_corrected(xenon_flash.raw_mspct, correction.method= MAYP11278_ylianttila.mthd, return.cps = TRUE)
getTimeUnit(xenon_flash.cps_spct)
## [1] "exposure"
plot(xenon_flash.cps_spct, range = c(315, NA))