1 Importing the output from TUV
As package ‘photobiologyInOut’ (versions >= 0.4.15) provides a function for the automatic import of the files returned by the Quick TUV Calculator, we show here how to produce such files, and link to an example showing how to import and plot them in R.
The TUV model is a well known model of atmospheric chemistry that also simulates the solar spectrum within the atmosphere or at ground level using a radiation transfer approach. The model is written in FORTRAN and available for local use. When one needs just to compute a few spectra with no special conditions, it is possible to use the on-line interface provided by the National Centre of Atmospheric Research (NCAR) under the name of Quick TUV Calculator.
This interface is easy to use once one understands the different options and required input data. (This video, authored by Pedro J. Aphalo is published with the express permission of Sasha Madronich, creator of the TUV model.)
This video is hosted at Vimeo.
2 Exploring sunlight with the Quick TUV model and R
I published an article describing how to import and plot the data from Quick TUV using R (Aphalo, 2018). A supplementary file to the article gives additional examples demonstrating how the solar spectrum is affected by variation in the atmospheric composition, solar elevation and observer elevation is available for reading and/or download. The R code used in embedded in the supplementary HTML file.
3 Diffuse and direct radiation under clear sky
[2024-08-27] Section added showing an example with multiple spectra, which also exemplifies how R package ‘photobiology’ and the R for photobiology suite makes it possible to work spectral data without having to code all the details of integration and interpolation. Package ‘ggspectra’ helps with the plotting of spectra.
The TUV model computes separately the diffuse and direct components of terrestrial solar spectral irradiance. I run 34 simulations varying only the sun elevation above the horizon. The aim was to show how the diffuse fraction depends both on sun elevation and on wavelengths under clear sky conditions. The files were manually saved from a web browser and ended as HTML files with leading and trailing lines added. These files can be read as is with ‘photobiologyInOut’ version >= 0.4.28.
I show here a full example, and start by removing all objects from the examples above.
We attach the packages that we will use.
We start by reading the 34 Quick TUV output files for the different solar elevation angles and split the components of total irradiance into separate spectral objects.
Code
# retrieve the file names
files <-
list.files(path = "QTUV-data", pattern = "tuv", full.names = TRUE)
# create an empty object to store the spectra
SZA.mspct <- source_mspct()
# read files one by one
for (f in files) {
# create short names from the file names
name <- gsub("tuv-default-|-4str.htm|-low.htm", "", basename(f))
print(name)
SZA.mspct[[name]] <- read_qtuv_txt(f, ozone.du = 300)
}
[1] "tuv-00.htm"
[1] "tuv-05.htm"
[1] "tuv-10.htm"
[1] "tuv-15.htm"
[1] "tuv-20.htm"
[1] "tuv-25.htm"
[1] "tuv-30.htm"
[1] "tuv-35.htm"
[1] "tuv-40.htm"
[1] "tuv-45.htm"
[1] "tuv-50.htm"
[1] "tuv-55.htm"
[1] "tuv-60.htm"
[1] "tuv-63.htm"
[1] "tuv-65.htm"
[1] "tuv-68.htm"
[1] "tuv-70.htm"
[1] "tuv-73.htm"
[1] "tuv-75.htm"
[1] "tuv-77.htm"
[1] "tuv-78.htm"
[1] "tuv-80.htm"
[1] "tuv-83.htm"
[1] "tuv-85.htm"
[1] "tuv-86.htm"
[1] "tuv-87.htm"
[1] "tuv-88.htm"
[1] "tuv-89.htm"
[1] "tuv-90.htm"
[1] "tuv-91.htm"
[1] "tuv-92.htm"
[1] "tuv-95.htm"
[1] "tuv-97.htm"
[1] "tuv-99.htm"
We can plot some spectra as examples.
In the next plot, the black line describes total spectral irradiance, with the sky blue area showing the diffuse component and the orange area the direct component. The diffuse or scattered light comes from the sky after being scattered in the atmosphere while direct radiation arrives directly from the sun without being scattered in the way.
Code
ggplot(SZA.mspct[["tuv-45.htm"]], aes(w.length, s.e.irrad)) +
geom_spct(fill = "skyblue", colour = "black") +
geom_spct(aes(y = s.e.irrad.dir), fill = "orange") +
scale_x_wl_continuous() +
scale_y_s.e.irrad_continuous()
Code
# each file contains three spectra that we separate here
SZA_tot.mspct <- source_mspct()
SZA_diff.mspct <- source_mspct()
SZA_dir.mspct <- source_mspct()
for (n in names(SZA.mspct)) {
SZA_tot.mspct[[n]] <-
SZA.mspct[[n]][ , c("w.length", "s.e.irrad")]
SZA_diff.mspct[[n]] <-
with(SZA.mspct[[n]],
source_spct(w.length = w.length,
s.e.irrad = s.e.irrad.diff.down))
SZA_dir.mspct[[n]] <-
with(SZA.mspct[[n]],
source_spct(w.length = w.length,
s.e.irrad = s.e.irrad.dir))
}
Compute irradiances and the diffuse fraction for different SZAs.
Code
wavebands <- c(list(PAR = PAR()), Plant_bands("sensory"))
wbands_tot.tb <-
q_irrad(SZA_tot.mspct, wavebands, scale.factor = 1e6)
wbands_dir.tb <-
q_irrad(SZA_dir.mspct, wavebands, scale.factor = 1e6)
wbands_diff.tb <-
q_irrad(SZA_diff.mspct, wavebands, scale.factor = 1e6)
wbands.tb <-
data.frame(
SZA = as.numeric(gsub("tuv-|\\.htm", "", wbands_tot.tb$spct.idx)),
SEA = 90 - as.numeric(gsub("tuv-|\\.htm", "", wbands_tot.tb$spct.idx))
)
for (col in setdiff(colnames(wbands_tot.tb), "spct.idx")) {
wbands.tb[[paste(col, "tot", sep = ".")]] <- wbands_tot.tb[[col]]
wbands.tb[[paste(col, "dir", sep = ".")]] <- wbands_dir.tb[[col]]
wbands.tb[[paste(col, "diff", sep = ".")]] <- wbands_diff.tb[[col]]
wbands.tb[[paste(col, "diff_fr", sep = ".")]] <-
wbands_diff.tb[[col]] / wbands_tot.tb[[col]]
}
# shorten column names
colnames(wbands.tb) <-
gsub("Q_]|Q_|\\.ISO|\\.CIE|\\.Sellaro|\\.Smith20", "", colnames(wbands.tb))
colnames(wbands.tb) # irradiances in umol m-2 s-1
[1] "SZA" "SEA" "PAR.tot" "PAR.dir"
[5] "PAR.diff" "PAR.diff_fr" "UVB.tot" "UVB.dir"
[9] "UVB.diff" "UVB.diff_fr" "UVA2.tot" "UVA2.dir"
[13] "UVA2.diff" "UVA2.diff_fr" "UVA1.tot" "UVA1.dir"
[17] "UVA1.diff" "UVA1.diff_fr" "Blue.tot" "Blue.dir"
[21] "Blue.diff" "Blue.diff_fr" "Green.tot" "Green.dir"
[25] "Green.diff" "Green.diff_fr" "Red.tot" "Red.dir"
[29] "Red.diff" "Red.diff_fr"
Plot PAR and UV-B normalized. We use solar elevation angle for the plot. The proportion of UV-B and to a lesser extent UV-A2 radiation is larger when the sun is higher above the horizon.
Code
ggplot(wbands.tb, aes(x = SEA)) +
geom_line(aes(y = PAR.tot / max(PAR.tot))) +
geom_line(aes(y = UVB.tot / max(UVB.tot)), colour = "purple") +
geom_line(aes(y = UVA2.tot / max(UVA2.tot)), colour = "violet") +
expand_limits(y = 0) +
labs(x = "Sun elevation angle (degrees)",
y = expression("Photon irradiance, "*italic(Q)~~("rel. units")))
In this example I was interested in the diffuse fraction. However, splines can be fitted to any of the components if desired.
Code
# We fit splines to the diffuse fraction obtaining functions
# that can be used to obtain by interpolation estimates for
# any solar elevation.
spl_funs.ls <- list()
diff_fr_cols <- grep("diff[_]fr", colnames(wbands.tb), value = TRUE)
for (col in diff_fr_cols) {
spl_funs.ls[[col]] <- splinefun(wbands.tb[["SEA"]], wbands.tb[[col]])
}
Another demonstration using PAR with the spline as a line and the simulated values as orange points.
Code
ggplot(wbands.tb, aes(SEA, PAR.diff_fr)) +
stat_function(fun = spl_funs.ls[["PAR.diff_fr"]], xlim = c(-10, 90), colour = "black") +
geom_point(na.rm = TRUE, colour = "orange") +
scale_x_continuous(name = "Solar elevation angle (degrees)",
breaks = c(-10, 0, 15, 30, 45, 60, 75, 90)) +
labs(y = expression("Diffuse fraction, "*Q[s] / Q[t]~~(""/1))) +
expand_limits(y = 0)
We have created a named list of function definitions.
Code
names(spl_funs.ls)
[1] "PAR.diff_fr" "UVB.diff_fr" "UVA2.diff_fr" "UVA1.diff_fr"
[5] "Blue.diff_fr" "Green.diff_fr" "Red.diff_fr"
A plot using the splines functions for single “colour” wavebands can now be easily created using the spline functions, which work like a normal function.
Code
ggplot(wbands.tb, aes(SEA)) +
stat_function(fun = spl_funs.ls[["UVB.diff_fr"]], xlim = c(-10, 90), colour = "black") +
stat_function(fun = spl_funs.ls[["UVA2.diff_fr"]], xlim = c(-10, 90), colour = "purple") +
stat_function(fun = spl_funs.ls[["UVA1.diff_fr"]], xlim = c(-10, 90), colour = "violet") +
stat_function(fun = spl_funs.ls[["Blue.diff_fr"]], xlim = c(-10, 90), colour = "blue") +
stat_function(fun = spl_funs.ls[["Green.diff_fr"]], xlim = c(-10, 90), colour = "green") +
stat_function(fun = spl_funs.ls[["Red.diff_fr"]], xlim = c(-10, 90), colour = "red") +
scale_x_continuous(name = "Solar elevation angle (degrees)",
breaks = c(-10, 0, 15, 30, 45, 60, 75, 90)) +
labs(y = expression("Diffuse fraction, "*Q[s] / Q[t]~~(""/1))) +
expand_limits(y = 0)
This plot shows clearly that solar radiation of longer wavelengths is less scattered and how scattering increases as the solar elevation decreases. Concurrently to the increased scattering, the equivalent length of the path through the atmosphere increases as shown below. The computaion used here is based on an emperirical formula that deviates from simple geometry.
Code
AM.df <- data.frame(SEA = 1:90,
AM = relative_AM(1:90))
ggplot(AM.df, aes(SEA, AM)) +
geom_line() +
scale_x_continuous(name = "Solar elevation angle (degrees)",
breaks = c(-10, 0, 15, 30, 45, 60, 75, 90)) +
labs(y = "Relative air mass, AM") +
expand_limits(y = 0)
4 Alternatives
R package ‘photobiologyInOut’ (>= 0.4.28-1) includes two functions that directly retrieve the spectral data files from the Quick TUV Calculator server and import them into R using the function used in the examples above. They support the same inputs as the web interface. Please, see the documentation.
The TUV model can be run on the local computer obtaining output files similar, but not identical, to those saved from the Quick TUV on-line calculator. Package ‘photobiologyInOut’ provides two functions to import these files. This would be preferable when the number of spectral simulations to do is large.
R package ‘foqat’ provides functions for calling on-line or locally the TUV model to produce batches of spectral simulations. The output is not directly compatible with the classes defined in package ‘photobiology’ to store spectra.