1 Introduction

The purpose of this R notebook is to demonstrate some of the recently added functionality of the R for Photobiology suite R for Photobiology suite of packages using a longer example than is possible to include in the documentation of the individual packages.

This file is an R notebook created using R markdown. It is self contained except for simulated spectral data which is imported from the output of the Quick TUV Calculator.

The code used to generate the HTML file is embedded in the HTML file itself as a script. The pretty-printed code chunks are interspersed with the text and figures, but hidden by default. There is a master button near the top-right corner to display all code, and individual buttons for each code chunk.

2 Use case example

This example was inspired by discussions at the Ultraviolet Photography discussion forum site on the solar spectrum.

Models are available for simulating the solar spectrum at ground level. The most frequently used ones are TUV and libradtran. Both models are radiative transfer models, i.e they use the extraterrestrial solar spectrum as input and based on the interaction of solar radiation with different layers of the atmosphere and with the Earth surface compute the solar spectrum at the surface or a other elevations within the atmosphere. Important characteristics of the atmosphere are the concentrations of different chemical species and the presence of different kinds of aerosols and clouds.

For intensive simulation exercises both models can be run locally. For smaller jobs, the on-line Quick TUV calculator is handy as I does not require any set up and parameters can be set through a form based interface. This on-line version gives access to only a subset of the functionality of TUV. For this exploration exercise, the settings are flexible enough.

The output of the results from the Quick TUV calculator, is simplified in its contents compared to that from the TUV program when run locally with output sent to a text file. In addition it may contain an HTML header and footer unless saved expressely as a text file.

library(readr)
library(photobiology)
library(photobiologyInOut)
library(photobiologySun)
library(photobiologyWavebands)
library(tibble)
library(tidyr)
library(ggplot2)
library(ggspectra)
library(dplyr)
library(lubridate)
# set the default theme for plots
theme_set(theme_bw(9))

In this article as published we plot spectral photon irradiances , but it can be easily rebuilt plotting spectral energy irradiances instead.

# Move the "#" to activate one of the two lines to reexpress all plots in the report
# energy_as_default()
photon_as_default()

The code below reads simulation resultss. It first searches in the current directory for all files with names starting with “tuv-” and ending with either “html”, “txt” or “text”, subsequently it reads them creating a collection of source spectra. This code is generally useful, and if a different naming convention is used, it is enough to change the search pattern in the first statement.

When we plot these spectra in later sections we index them by name. The chunk above names the members of the collection based on file names. Dashes (“-”) are replaced by underscores ("_") to make names that do not require quoting when used. If additional special characters are present in file names, you may want to use R function make.names() to make syntactically valid names out of the file names.

3 Extraterrestrial

As a reference we have in the first plot the extraterrestrial, or zero air mass (AM0) solar spectrum according to Wehrli as used by the World Meteorological Organization (WMO). All peaks and valleys visible in this plot are related to emission and absorption phenomena taking place outside the Earth atmosphere.

autoplot(sun_reference.mspct$WMO.Wehrli.AM0, 
     range = c(NA, 780), 
     annotations = c("-", "peaks"))

4 Standard terrestrial

The ASTM has published a standard spectrum for one and half air mass (AM1.5). It is based on simulations and attempts to represent the long-term average conditions in the USA. In this figure it is possible to clearly see some deep valleys caused be absorption by atmospheric gases.

autoplot(sun_reference.mspct$ASTM.G173.global, 
     range = c(NA, 780), 
     annotations = c("-", "peaks"))

5 Solar elevation

Instead of thinking in terms of seasons and latitudes, considering the solar elevation angle gives us more generally applicable simulations. Spectral photon irradiance decreases with decreasing sun elevation angle, but differences in the shape of the spectra are difficult to fathom when the spectra are plotted in absolute units.

autoplot(sun_tuv.mspct[c("tuv_azimuth_00_O3_300",
                     "tuv_azimuth_15_O3_300",
                     "tuv_azimuth_30_O3_300",
                     "tuv_azimuth_45_O3_300",
                     "tuv_azimuth_60_O3_300",
                     "tuv_azimuth_75_O3_300",
                     "tuv_azimuth_85_O3_300",
                     "tuv_azimuth_88_O3_300",
                     "tuv_azimuth_90_O3_300"
                     )], 
     annotations = c("-", "peaks")) +
  scale_linetype(labels = 90 - c(0, 15, 30, 45, 60, 75, 85, 88, 90),
                 name = "Sun\nelevation\n(degrees)")

Even a detail of the UV region remains difficult to “read”.

autoplot(sun_tuv.mspct[c("tuv_azimuth_00_O3_300",
                     "tuv_azimuth_15_O3_300",
                     "tuv_azimuth_30_O3_300",
                     "tuv_azimuth_45_O3_300",
                     "tuv_azimuth_60_O3_300",
                     "tuv_azimuth_75_O3_300",
                     "tuv_azimuth_85_O3_300",
                     "tuv_azimuth_88_O3_300",
                     "tuv_azimuth_90_O3_300"
                     )], 
     range = c(280, 380),
     annotations = c("-", "peaks")) +
  scale_linetype(labels = 90 - c(0, 15, 30, 45, 60, 75, 85, 88, 90),
                 name = "Sun\nelevation\n(degrees)")

If we divide all spectra by that for the sun at the zenith (elevation angle of 90 degrees) we obtain an easier to interpret plot, where we can see that the proportion of UVB photons decreases very quickly when the sun is lower than at the zenith.

convolve_each(
  sun_tuv.mspct[c("tuv_azimuth_00_O3_300",
                  "tuv_azimuth_15_O3_300",
                  "tuv_azimuth_30_O3_300",
                  "tuv_azimuth_45_O3_300",
                  "tuv_azimuth_60_O3_300",
                  "tuv_azimuth_75_O3_300",
                  "tuv_azimuth_85_O3_300",
                  "tuv_azimuth_88_O3_300",
                  "tuv_azimuth_90_O3_300"
  )], sun_tuv.mspct[["tuv_azimuth_00_O3_300"]], oper = `/`) %>%
  msmsply(setScaled, scaled = TRUE) %>%
plot(annotations = c("-", "peaks")) +
  scale_linetype(labels = 90 - c(0, 15, 30, 45, 60, 75, 85, 88, 90),
                 name = "Sun\nelevation\n(degrees)")

Another interesting question is the contribution of diffuse radiation to total (diffuse + direct) spectral photon irradiance and its dependence on elevation angle and wavelength. Clearly the fraction of diffuse radiation increases as the sun approaches the horizon, but in addition the fraction of diffuse radiation is much larger in the UV than in the red and NIR regions. This is the reason why contrast between shadows and full sun is much less in UV than in VIS. As NIR is even less scattered distant objects are clearer in NIR than in VIS bands.

ggplot(sun_tuv.mspct[c("tuv_azimuth_00_O3_300",
                     "tuv_azimuth_15_O3_300",
                     "tuv_azimuth_30_O3_300",
                     "tuv_azimuth_45_O3_300",
                     "tuv_azimuth_60_O3_300",
                     "tuv_azimuth_75_O3_300",
                     "tuv_azimuth_85_O3_300",
                     "tuv_azimuth_88_O3_300",
                     "tuv_azimuth_90_O3_300"
                     )], 
     aes(w.length, s.e.irrad.diff.down/s.e.irrad, linetype = spct.idx)) +
  geom_line(na.rm = TRUE) +
  scale_linetype(labels = 90 - c(0, 15, 30, 45, 60, 75, 85, 88, 90),
                 name = "Sun\nelevation\n(degrees)") +
  scale_x_wl_continuous() +
  scale_colour_identity() +
  expand_limits(y = 0) +
  labs(y = "Diffuse radiation, fraction of global radiation (/1)") +
  ggtitle("Diffuse radiation under clear sky")

6 Ozone column

Being interested in UV radiation, the obvious next question to tackle is ozone column depth. In the next simulation, for the Equator and using very extreme values for the ozone column we see a change in non-weighted UV-B radiation is difficult to see in the whole spectrum plot.

plot(sun_tuv.mspct[c("tuv_azimuth_00_O3_400", 
                     "tuv_azimuth_00_O3_300", 
                     "tuv_azimuth_00_O3_200")], 
     annotations = c("-", "peaks")) +
  scale_linetype(labels = c(400, 300, 200),
                 name = "Ozone\ncolum\n(d.u.)")

Looking at a detail of the UV region we can see that the effect in the UV-A is barely noticeable, and that in the UV-B it is larger. This effect looks still surprisingly small given the huge range of ozone column values used for the simulation. The integrated values for each band in the UV are displayed on the bars.

plot(sun_tuv.mspct[c("tuv_azimuth_00_O3_400", 
                     "tuv_azimuth_00_O3_300", 
                     "tuv_azimuth_00_O3_200")],
     range = c(NA, 380), 
     idfactor =  NA, 
     annotations = c("-", "peaks")) + 
  facet_wrap(~spct.idx)

The figure above considers only the total number of photons, ignoring that photons at shorter wavelength can have a stronger effect on organisms. If we replot the same data weighted with erythremal spectrum (for the reddening of human skin) the expected effect of ozone depletion on human skin becomes easier to grasp. In this case we use energy units as this is customary when calculating biologically effective irradiances.

convolve_each(
  sun_tuv.mspct[c("tuv_azimuth_00_O3_400", 
                  "tuv_azimuth_00_O3_300", 
                  "tuv_azimuth_00_O3_200")], CIE()) %>%
  plot(range = c(NA, 380), 
       idfactor =  NA, 
       annotations = c("-", "peaks"),
       unit.out = "energy") + 
  facet_wrap(~spct.idx)

7 Elevation above sea level

An increase in the elevation where the observer is situated reduces the air mass in the pass of the solar radiation reaching this observer. In the next simulations we assume that the observer is at ground level but that the ground surface is at different elevations above sea level. The increase in fine structure is because I did this simulation using a wavelength step of 0.2 nm. The most frequently used spectroradiometers have lower wavelength resolution. In reality the solar spectrum has even finer structure.

plot(sun_tuv.mspct[c("tuv_azimuth_00_O3_300_hires",
                     "tuv_azimuth_00_O3_300_hires_2km",
                     "tuv_azimuth_00_O3_300_hires_4km")],
     idfactor =  NA, 
     annotations = c("-", "peaks")) + 
  facet_wrap(~spct.idx)

convolve_each(
  sun_tuv.mspct[c("tuv_azimuth_00_O3_300_hires",
                  "tuv_azimuth_00_O3_300_hires_2km",
                  "tuv_azimuth_00_O3_300_hires_4km")], CIE()) %>%
  plot(range = c(NA, 380), 
       idfactor =  NA, 
       annotations = c("-", "peaks")) + 
  facet_wrap(~spct.idx)

8 Cloud depth

Clouds affect both the irradiance at ground level through absorption and reflection, and increase the proportion of diffuse radiation in global radiation as a result of scattering. Precisely describing clouds, and simulating their effect on spectral irradiance is difficult. Here we use four arbitrary values for optical cloud depth.

plot(sun_tuv.mspct[c("tuv_azimuth_45_O3_300_clouds00",
                     "tuv_azimuth_45_O3_300_clouds25",
                     "tuv_azimuth_45_O3_300_clouds50",
                     "tuv_azimuth_45_O3_300_clouds100")],
     idfactor =  NA, 
     annotations = c("-", "peaks")) + 
  facet_wrap(~spct.idx, nrow = 1L)

convolve_each(
  sun_tuv.mspct[c("tuv_azimuth_45_O3_300_clouds00",
                  "tuv_azimuth_45_O3_300_clouds25",
                  "tuv_azimuth_45_O3_300_clouds50",
                  "tuv_azimuth_45_O3_300_clouds100")], CIE()) %>%
  plot(range = c(NA, 380), 
       idfactor =  NA, 
       annotations = c("-", "peaks")) + 
  facet_wrap(~spct.idx, nrow = 1L)

9 Appendix: Solar elevation angles

Here we present plots of solar elevation through the course of the day, at the solstices at six different latitudes.

9.1 Summer solstice

sun_angles(time = ymd_hms("2018-06-21 00:00:00", tz = "UTC") + minutes(0:(23 * 60)), 
           geocode = data.frame(lon = 0, lat = c(0, 16.5, 33, 49.5, 66, 82.5) ) ) %>%
  ggplot(aes(solartime, elevation, colour = factor(latitude))) +
  annotate(geom = "rect", 
           xmin = -Inf, xmax = Inf, ymin = -90, ymax = 0, 
           colour = "black", alpha = 0.1, size  = 0,
           linetype = "dashed") +
  geom_hline(yintercept = 90, linetype = "dashed") +
  geom_line() +
  expand_limits(y = c(-90,90)) +
  scale_x_continuous(breaks = c(0,6,12,18,24)) +
  scale_y_continuous(breaks = c(-90, -45, 0, 45, 90)) +
  labs(y = "Sun elevation angle (degrees above the horizon)",
       x = "Time of day, local solar time (h)",
       colour = "Latitude\n(degrees)")