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 or in the Bulletin article to which this file is a supplement.
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, shown by default. There is a master button near the top-right corner to hide all R code, and individual buttons for each code chunk.
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 results. 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.
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.
spectra7elevations <- c("tuv_zenith_00_O3_300",
"tuv_zenith_15_O3_300",
"tuv_zenith_30_O3_300",
"tuv_zenith_45_O3_300",
"tuv_zenith_60_O3_300",
"tuv_zenith_75_O3_300",
"tuv_zenith_90_O3_300")
the7elevations <- 90 - c(0, 15, 30, 45, 60, 75, 90)
autoplot(sun_tuv.mspct[spectra7elevations],
annotations = c("-", "peaks")) +
scale_linetype(labels = the7elevations,
name = "Sun\nelevation\n(degrees)")
Even a detail of the UV region remains difficult to “read”.
autoplot(sun_tuv.mspct[spectra7elevations],
range = c(280, 380),
annotations = c("-", "peaks")) +
scale_linetype(labels = the7elevations,
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[spectra7elevations],
sun_tuv.mspct[["tuv_zenith_00_O3_300"]],
oper = `/`) %>%
msmsply(setScaled, scaled = TRUE) %>%
plot(annotations = c("-", "peaks")) +
scale_linetype(labels = the7elevations,
name = "Sun\nelevation\n(degrees)")
q_ratio(sun_tuv.mspct[spectra7elevations],
list(UVB(), UVA(), Blue("Sellaro")), PAR())
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[spectra7elevations],
aes(w.length, s.e.irrad.diff.down/s.e.irrad, linetype = spct.idx)) +
geom_line(na.rm = TRUE) +
scale_linetype(labels = the7elevations,
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")
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_zenith_00_O3_400",
"tuv_zenith_00_O3_300",
"tuv_zenith_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_zenith_00_O3_400",
"tuv_zenith_00_O3_300",
"tuv_zenith_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_zenith_00_O3_400",
"tuv_zenith_00_O3_300",
"tuv_zenith_00_O3_200")], CIE()) %>%
plot(range = c(NA, 380),
idfactor = NA,
annotations = c("-", "peaks"),
unit.out = "energy") +
facet_wrap(~spct.idx)
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_zenith_45_O3_300_clouds00",
"tuv_zenith_45_O3_300_clouds25",
"tuv_zenith_45_O3_300_clouds50",
"tuv_zenith_45_O3_300_clouds100")],
idfactor = NA,
annotations = c("-", "peaks")) +
facet_wrap(~spct.idx, nrow = 1L)
convolve_each(
sun_tuv.mspct[c("tuv_zenith_45_O3_300_clouds00",
"tuv_zenith_45_O3_300_clouds25",
"tuv_zenith_45_O3_300_clouds50",
"tuv_zenith_45_O3_300_clouds100")], CIE()) %>%
plot(range = c(NA, 380),
idfactor = NA,
annotations = c("-", "peaks")) +
facet_wrap(~spct.idx, nrow = 1L)
Above we presented plots of solar elevation through the course of the day, at the solstices at six different latitudes. In this appendix we show how to compute and plot the angles defining the position of the sun.
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)")
sun_angles(time = ymd_hms("2018-03-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)")
sun_angles(time = ymd_hms("2018-12-21 01: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)")
The air mass that the solar beam travels through depends on the elevation angle, the higher the sun is in the sky, the shorter the path. When the sun is at the zenith the beam crosses the atmosphere perpendicularly and consequently travels through one depth of the atmosphere (AM1).
One possible approximation is geometric, based on the diameter of the Earth and the depth of the atmosphere. Empirical formulas have been developed that better describe the effective length of the path. The different approaches result is somehow different values at low sun elevations.
tibble(elevation.angle = 1:90, AM = relative_AM(1:90)) %>%
ggplot(aes(elevation.angle, AM)) +
geom_line() +
geom_hline(yintercept = 1, linetype = "dashed") +
scale_x_continuous(breaks = c(0, 10, 20, 45, 90), minor_breaks = numeric()) +
scale_y_continuous(breaks = c(1, 3, 5, 10, 15, 20, 25), minor_breaks = numeric()) +
labs(y = "Air Mass, AM (reltive units)",
x = "Sun elevation angle (degrees above the horizon)")
Here we automatically list the R version and package versions used to create the file you are reading.
- Session info -------------------------------------------------------------------------
- Packages -----------------------------------------------------------------------------
package * version date lib source
assertthat 0.2.0 2017-04-11 [1] CRAN (R 3.5.0)
backports 1.1.2 2017-12-13 [1] CRAN (R 3.5.0)
bindr 0.1.1 2018-03-13 [1] CRAN (R 3.5.0)
bindrcpp 0.2.2 2018-03-29 [1] CRAN (R 3.5.0)
callr 3.1.0 2018-12-10 [1] CRAN (R 3.5.1)
cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.1)
colorspace 1.3-2 2016-12-14 [1] CRAN (R 3.5.0)
colorSpec 0.7-5 2018-11-19 [1] CRAN (R 3.5.1)
crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.0)
desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.0)
devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.1)
digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.1)
dplyr * 0.7.8 2018-11-10 [1] CRAN (R 3.5.1)
fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.1)
ggplot2 * 3.1.0 2018-10-25 [1] CRAN (R 3.5.1)
ggrepel 0.8.0 2018-05-09 [1] CRAN (R 3.5.1)
ggspectra * 0.3.1 2018-10-29 [1] local
glue 1.3.0 2018-07-17 [1] CRAN (R 3.5.1)
gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.0)
hms 0.4.2 2018-03-10 [1] CRAN (R 3.5.0)
knitr 1.21 2018-12-10 [1] CRAN (R 3.5.1)
labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.0)
lubridate * 1.7.4 2018-04-11 [1] CRAN (R 3.5.0)
magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.0)
MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.1)
memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.0)
minpack.lm 1.2-1 2016-11-20 [1] CRAN (R 3.5.0)
munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.0)
photobiology * 0.9.25 2018-12-09 [1] local
photobiologyInOut * 0.4.19 2018-12-13 [1] local
photobiologySun * 0.4.0.9000 2018-06-10 [1] local
photobiologyWavebands * 0.4.2-1.9001 2018-10-18 [1] local
pillar 1.3.0.9001 2018-11-28 [1] Github (r-lib/pillar@c5bf622)
pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.1)
pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.1)
plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.0)
prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.0)
processx 3.2.1 2018-12-05 [1] CRAN (R 3.5.1)
ps 1.2.1 2018-11-06 [1] CRAN (R 3.5.1)
purrr 0.2.5 2018-05-29 [1] CRAN (R 3.5.0)
R6 2.3.0 2018-10-04 [1] CRAN (R 3.5.1)
Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.1)
readr * 1.3.0 2018-12-11 [1] CRAN (R 3.5.1)
remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.1)
rlang 0.3.0.1 2018-10-25 [1] CRAN (R 3.5.1)
rootSolve 1.7 2016-12-06 [1] CRAN (R 3.5.0)
rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.0)
rstudioapi 0.8 2018-10-02 [1] CRAN (R 3.5.1)
scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.1)
stringi 1.2.4 2018-07-20 [1] CRAN (R 3.5.1)
stringr 1.3.1 2018-05-10 [1] CRAN (R 3.5.0)
testthat 2.0.1 2018-10-13 [1] CRAN (R 3.5.1)
tibble * 1.4.99.9006 2018-12-03 [1] Github (tidyverse/tibble@ddc9110)
tidyr * 0.8.2 2018-10-28 [1] CRAN (R 3.5.1)
tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.1)
withr 2.1.2 2018-08-27 [1] Github (jimhester/withr@8b9cee2)
xfun 0.4 2018-10-23 [1] CRAN (R 3.5.1)
yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1)
[1] C:/Users/aphalo/Documents/R/win-library/3.5
[2] C:/Program Files/R/R-3.5.1patched/library
Copyright (c) 2018 Pedro J. Aphalo.