Advanced Examples
Package ‘photobiologySun’ 0.5.0.9000
Pedro J. Aphalo
2024-09-17
Source:vignettes/articles/advanced-examples.Rmd
advanced-examples.Rmd
library(photobiology)
library(photobiologyWavebands)
library(photobiologyInOut)
library(photobiologySun)
library(lubridate)
library(ggspectra)
eval_plots <- requireNamespace("ggspectra", quietly = TRUE)
if (eval_plots) {
library(ggspectra)
theme_set(theme_bw())
}Introduction
This article is included in the on-line documentation, but not in the package distribution, as doing its rendering a large file is downloaded. It also requires additional R packages to be installed.
High frequency irradiance data
The irradiance data in object four_days_1min.data are
only for four days out of a time-series of data started in 2017 and
still continuing (Aphalo, 2023). The whole data set can be downloaded at
doi:10.17605/OSF.IO/E4VAU
or https://osf.io/e4vau/. Data
logged at 1 min intervals or at 1 h intervals are available. Objects
ppfd_LICOR.data, ppfd_BF.data and
irrad_Kipp.data contain earlier data from the same weather
station for a period of 17 days.
At the moment I do not know how to retrieve the download URL of the individual data files programmatically. It can be found through the web interface, or the files can be manually downloaded through this same interface. The example in this section does the file download directly in R.
Files are large, so in many cases it will be necessary to increase
the value of the timeout R option temporarily to several
minutes from its default of 60 s.
When the download time-outs before completion, a smaller and corrupted file is saved to disk. The code below downloads the file only once. In the case of a truncated or corrupted file a new download requires the local file to be deleted.
old.timeout <- options(timeout = 180L)
if (!file.exists("Viikki_1min_2023.rda")) {
download.file("https://osf.io/download/3rbd9/",
destfile = "Viikki_1min_2023.rda",
cacheOK = TRUE,
mode = "wb")
}
options(old.timeout)We can now load the file into the R workspace.
load("Viikki_1min_2023.rda")
class(Viikki_1min_2023_Version_1.0.tb)
#> [1] "tbl_df" "tbl" "data.frame"
ncol(Viikki_1min_2023_Version_1.0.tb)
#> [1] 66
nrow(Viikki_1min_2023_Version_1.0.tb)
#> [1] 525567The file contains data for one year.
range(Viikki_1min_2023_Version_1.0.tb$time)
#> [1] "2023-01-01 00:00:00 UTC" "2023-12-31 23:59:00 UTC"The file contains in addition to solar radiation, data for other weather variables, both measured and derived by calculation.
cat(colnames(Viikki_1min_2023_Version_1.0.tb), sep = "\n")
#> time
#> day_of_year
#> month_of_year
#> month_name
#> calendar_year
#> time_of_day
#> solar_time
#> sun_elevation
#> sun_azimuth
#> PAR_umol_LI
#> PAR_umol_CS
#> PAR_umol_BF
#> PAR_umol
#> PAR_diff_fr
#> global_watt
#> red_umol
#> far_red_umol
#> blue_umol
#> blue_sellaro_umol
#> UVA_umol
#> UVA1_umol
#> UVA2_umol
#> UVB_umol
#> blue_red
#> blue_red_sq
#> UVA_PAR
#> UVA_PAR_sq
#> UVA1_PAR
#> UVA1_PAR_sq
#> UVA2_PAR
#> UVA2_PAR_sq
#> UVB_PAR
#> UVB_PAR_sq
#> red_far_red
#> wind_speed
#> wind_direction
#> air_temp_C
#> air_temp_min_C
#> air_temp_max_C
#> air_vp
#> air_RH
#> air_DP
#> air_pressure
#> rain_mm_min
#> surf_temp_C
#> surf_temp_sensor_delta_C
#> was_sunny
#> solar_time_h
#> was_day
#> PAR_diff_fr_rel
#> logged_air_temp_C
#> air_temp_run_median
#> temp_surf2air_C
#> R_0
#> R_rel
#> Rn_sw_ref
#> Rn_ref
#> ET_ref_FAO56
#> ET_ref_short
#> ET_ref_tall
#> air_vpd
#> SupplyVoltage_Max
#> ReferenceVoltage_Min
#> ReferenceVoltage_Max
#> BattV_Min
#> BattV_Max
cat(comment(Viikki_1min_2023_Version_1.0.tb))
#> Data acquired or computed by Pedro J. Aphalo
#> Data are provided as is, with no guaranee of suitability for any purpose
#> If used in publications cite based on DOI
#> Faculty of Biological and Environmental Sciences, University of Helsinki, Finland
#> License: CC BY-NC-SA Attribution-NonCommercial-ShareAlike 4.0 International
#> See: https://creativecommons.org/licenses/by-nc-sa/4.0/
#> Station information available at https://viikki-stn.r4photobiology.info/
#> Period: 2023-01-01 to 2023-12-31 23:59:00
#> Filename: Viikki-1min-2023-Version-1.0.rda
print(object.size(Viikki_1min_2023_Version_1.0.tb),
standard = "SI", units = "MB")
#> 271.2 MBGiven the size of the object, it is best to subset rows and select
columns that are of interest before plotting or computing summaries from
the data set. This is specially important for ‘ggplot2’ as a copy of the
argument passed to data is added to the ggplot
("gg") object. Here we plot data for a single day, against
local solar time.
subset(Viikki_1min_2023_Version_1.0.tb,
subset = time >= ymd_hms("2023-6-21 00:00:00") &
time < ymd_hms("2023-6-22 00:00:00"),
select = c(solar_time, PAR_umol)) |>
ggplot(aes(solar_time, PAR_umol)) +
geom_line()
subset(Viikki_1min_2023_Version_1.0.tb,
subset = time >= ymd_hms("2023-6-21 00:00:00") &
time < ymd_hms("2023-6-22 00:00:00"),
select = c(solar_time, global_watt)) -> one_day.tb
one_day.tb[["global_watt"]] <-
ifelse(one_day.tb[["global_watt"]] < 0, 0, one_day.tb[["global_watt"]])
daily_MJ_on_site <- sum(one_day.tb[["global_watt"]]) * 60 * 1e-6 # MJ m-2 d-1
daily_MJ_on_site
#> [1] 26.83269
subset(Viikki_1min_2023_Version_1.0.tb,
subset = time >= ymd_hms("2023-6-21 00:00:00") &
time < ymd_hms("2023-6-22 00:00:00"),
select = c(solar_time, global_watt)) |>
ggplot(aes(solar_time, global_watt)) +
geom_line() +
expand_limits(y = 900) +
labs(y = expression("Global radiation "*(W~m^{-2})),
x = "Solar time (h)")
If we do no longer need the downloaded file, it can be deleted running the code below (here not run).
unlink("Viikki_1min_2023.rda")Other sources of irradiance data include: NASA Power (see R packages ‘nasapower’ and ‘photobiologyInOut’), with worldwide coverage at ground level; SORCE, with daily spectral data for solar radiation at the top of the atmosphere; EUMETSAT CM SAF includes surface radiation products for climatology but no real-time data; EUMETSAT CAMS provides global radiation and its components with a delay of 2 d with coverage of Europe and nearby regions (see R package ‘camsRad’), data accessible through SoDa.
NASA Power data
Data available through NASA Power are mainly intended for estimates of renewable energy generation, but also include some meteorological data. The API supports arbitrary geographic coordinates as input.
library(nasapower)
pars.available <- query_parameters(community = "ag", temporal_api = "hourly")
str(pars.available, max.level = 1)
#> List of 52
#> $ PRECSNOLAND :List of 9
#> $ PRECTOTCORR :List of 9
#> $ PS :List of 9
#> $ QV10M :List of 9
#> $ QV2M :List of 9
#> $ RH2M :List of 9
#> $ SNODP :List of 9
#> $ T2M :List of 9
#> $ TS :List of 9
#> $ U10M :List of 9
#> $ U2M :List of 9
#> $ U50M :List of 9
#> $ V10M :List of 9
#> $ V2M :List of 9
#> $ V50M :List of 9
#> $ PSC :List of 9
#> $ T2MDEW :List of 9
#> $ T2MWET :List of 9
#> $ WD10M :List of 9
#> $ WD2M :List of 9
#> $ WD50M :List of 9
#> $ WS10M :List of 9
#> $ WS2M :List of 9
#> $ WS50M :List of 9
#> $ WSC :List of 9
#> $ ALLSKY_SFC_LW_DWN :List of 9
#> $ ALLSKY_SFC_SW_DIFF :List of 9
#> $ ALLSKY_SFC_SW_DWN :List of 9
#> $ ALLSKY_SFC_UV_INDEX:List of 9
#> $ ALLSKY_SFC_UVA :List of 9
#> $ ALLSKY_SFC_UVB :List of 9
#> $ AOD_55 :List of 9
#> $ AOD_84 :List of 9
#> $ CLOUD_AMT :List of 9
#> $ CLOUD_OD :List of 9
#> $ CLRSKY_SFC_LW_DWN :List of 9
#> $ CLRSKY_SFC_SW_DIFF :List of 9
#> $ CLRSKY_SFC_SW_DWN :List of 9
#> $ PW :List of 9
#> $ SZA :List of 9
#> $ TOA_SW_DWN :List of 9
#> $ ALLSKY_KT :List of 9
#> $ ALLSKY_NKT :List of 9
#> $ ALLSKY_SFC_PAR_TOT :List of 9
#> $ ALLSKY_SFC_SW_DNI :List of 9
#> $ ALLSKY_SRF_ALB :List of 9
#> $ CLRSKY_KT :List of 9
#> $ CLRSKY_NKT :List of 9
#> $ CLRSKY_SFC_PAR_TOT :List of 9
#> $ CLRSKY_SFC_SW_DNI :List of 9
#> $ CLRSKY_SRF_ALB :List of 9
#> $ TOA_SW_DNI :List of 9
kumpula.global_watt.tb <-
get_power(community = "ag",
pars = c("ALLSKY_SFC_SW_DWN", "CLRSKY_SFC_SW_DWN", "ALLSKY_SFC_PAR_TOT", "CLRSKY_SFC_PAR_TOT"),
temporal_api = "hourly",
lonlat = c(25.01673, 60.2253), # Viikki, Helsinki, Finland
dates = "2023-6-21",
time_standard = "LST"
)
colnames(kumpula.global_watt.tb)
#> [1] "LON" "LAT" "YEAR"
#> [4] "MO" "DY" "HR"
#> [7] "ALLSKY_SFC_SW_DWN" "CLRSKY_SFC_SW_DWN" "ALLSKY_SFC_PAR_TOT"
#> [10] "CLRSKY_SFC_PAR_TOT"
str(pars.available$ALLSKY_SFC_SW_DWN)
#> List of 9
#> $ type : chr "RADIATION"
#> $ temporal : chr "HOURLY"
#> $ source : chr "CERES"
#> $ community : chr "AG"
#> $ calculated: logi FALSE
#> $ inputs : NULL
#> $ units : chr "MJ/hr"
#> $ name : chr "All Sky Surface Shortwave Downward Irradiance"
#> $ definition: chr "The total solar irradiance incident (direct plus diffuse) on a horizontal plane at the surface of the earth und"| __truncated__
kumpula.global_watt.tb[["global_watt"]] <-
kumpula.global_watt.tb[["ALLSKY_SFC_SW_DWN"]] * 1e6 / 3.6e3
daily_MJ_NASA_Power <- sum(kumpula.global_watt.tb[["global_watt"]]) * 3600 * 1e-6 # MJ m-2 d-1
daily_MJ_NASA_Power
#> [1] 24.53
daily_MJ_NASA_Power / daily_MJ_on_site
#> [1] 0.9141835
ggplot(kumpula.global_watt.tb, aes(HR, global_watt)) +
expand_limits(y = 900) +
geom_line() +
labs(y = expression("Global radiation "*(W~m^{-2})),
x = "Solar time (h)")
str(pars.available$ALLSKY_SFC_PAR_TOT)
#> List of 9
#> $ type : chr "RADIATION"
#> $ temporal : chr "HOURLY"
#> $ source : chr "POWER"
#> $ community : chr "AG"
#> $ calculated: logi FALSE
#> $ inputs : NULL
#> $ units : chr "W/m^2"
#> $ name : chr "All Sky Surface PAR Total"
#> $ definition: chr "The total Photosynthetically Active Radiation (PAR) incident on a horizontal plane at the surface of the earth "| __truncated__
ggplot(kumpula.global_watt.tb, aes(HR,
ALLSKY_SFC_PAR_TOT)) +
geom_line() +
labs(y = expression("PhR "*(W~m^{-2})),
x = "Solar time (h)")
EUMETSAT CAMS radiation service
Data from the EUMETSAT CAMS radiation service can be downloaded interactively or using API calls. The global radiation data available through this service is rather similar to that in NASA Power. However, one important difference is that data are natively at 1 min intervals in the time series, although also available for download integrated over longer time periods. Data are available both in .CSV (using English conventions for separators and decimal mark) and in netDCF format. To be able to download data, a free account should be created. For this example I downloaded the data interactively as a netCDF file. The maximum number of queries per day and user is enforced to 100. A single query can include many grid points and times. For this example I downloaded data for the same day as in the examples above, with a 1 min time resolution.
Viikki_CAMS.tnc <- tidync("Viikki-CAMS-solar-radiation-1min.nc")
Viikki_CAMS.tnc
#>
#> Data Source (1): Viikki-CAMS-solar-radiation-1min.nc ...
#>
#> Grids (5) <dimension family> : <associated variables>
#>
#> [1] D2,D1,D0,D3 : ut_jd, ut_year, ut_month, ut_day, ut_hour, ut_minute, ut_second, tst_jd, tst_year, tst_month, tst_day, tst_hour, tst_minute, tst_second, G0, CLEAR_SKY_GHI, CLEAR_SKY_BHI, CLEAR_SKY_DHI, CLEAR_SKY_BNI, rely, GHI, BHI, DHI, BNI **ACTIVE GRID** ( 1440 values per variable)
#> [2] D0 : altitude
#> [3] D1 : latitude
#> [4] D2 : longitude
#> [5] D3 : time
#>
#> Dimensions 4 (all active):
#>
#> dim name length min max start count dmin dmax unlim coord_dim
#> <chr> <chr> <dbl> <dbl> <dbl> <int> <int> <dbl> <dbl> <lgl> <lgl>
#> 1 D0 altitude 1 0 0 1 1 0 0 FALSE TRUE
#> 2 D1 latitude 1 6.02e1 6.02e1 1 1 6.02e1 6.02e1 FALSE TRUE
#> 3 D2 longitude 1 2.50e1 2.50e1 1 1 2.50e1 2.50e1 FALSE TRUE
#> 4 D3 time 1440 1.69e9 1.69e9 1 1440 1.69e9 1.69e9 FALSE TRUEFunction hyper_dims() returns a description of the grid
for which observations are available.
hyper_dims(Viikki_CAMS.tnc)
#> # A tibble: 4 × 7
#> name length start count id unlim coord_dim
#> <chr> <dbl> <int> <int> <int> <lgl> <lgl>
#> 1 longitude 1 1 1 2 FALSE TRUE
#> 2 latitude 1 1 1 1 FALSE TRUE
#> 3 altitude 1 1 1 0 FALSE TRUE
#> 4 time 1440 1 1440 3 FALSE TRUEFunction hyper_vars() returns a description of the
observations or variables available at each grid point.
hyper_vars(Viikki_CAMS.tnc)
#> # A tibble: 24 × 6
#> id name type ndims natts dim_coord
#> <int> <chr> <chr> <int> <int> <lgl>
#> 1 4 ut_jd NC_DOUBLE 4 2 FALSE
#> 2 5 ut_year NC_SHORT 4 2 FALSE
#> 3 6 ut_month NC_BYTE 4 2 FALSE
#> 4 7 ut_day NC_BYTE 4 2 FALSE
#> 5 8 ut_hour NC_BYTE 4 2 FALSE
#> 6 9 ut_minute NC_BYTE 4 2 FALSE
#> 7 10 ut_second NC_DOUBLE 4 2 FALSE
#> 8 11 tst_jd NC_DOUBLE 4 2 FALSE
#> 9 12 tst_year NC_SHORT 4 2 FALSE
#> 10 13 tst_month NC_BYTE 4 2 FALSE
#> # ℹ 14 more rowsFunction hyper_tibble() extracts a subset of the data
into a tibble in long (or “tidy”) format. The selection of the grid
point can be done in the same operation. A pipe is used to add the
decoded dates, using the pipe operator (|>) and methods from ‘dplyr’.
The decoding of dates is done using functions from package
‘lubridate’.
Viikki_CAMS.tb <- hyper_tibble(Viikki_CAMS.tnc)
Viikki_CAMS.tb[["solar.time"]] <-
with(Viikki_CAMS.tb, tst_hour + tst_minute / 60 + tst_second / 3600)
Viikki_CAMS.tb[["global.watt"]] <-
with(Viikki_CAMS.tb, GHI * 60)We can calculate the daily global radiation and compare it to that measured on site. The difference of less than 3% is well within the calibration tolerance of the pyranometer used.
daily_MJ_CAMS <- sum(Viikki_CAMS.tb[["global.watt"]]) * 60 * 1e-6 # MJ m-2 d-1
daily_MJ_CAMS
#> [1] 27.57183
daily_MJ_CAMS / daily_MJ_on_site
#> [1] 1.027546The time-course matches closely that measured on-site, with clouds near 5:00 local solar time. This is an incredibly good match for a satellite-based estimate to a ISO 9060 spectrally flat Class C pyranometer, designed for routine measurements.
Viikki_CAMS.tb |>
subset(select = c(solar.time, global.watt)) |>
ggplot(aes(solar.time, global.watt)) +
geom_line() +
expand_limits(y = 900) +
labs(y = expression("Global radiation "*(W~m^{-2})),
x = "Solar time (h)")
Interpolation functions for diffuse fraction
In this section we will use the example spectral irradiance data in
sun_elevation.spct obtained as a simulation with the TUV
model under clear sky conditions. We will compute the diffuse fraction
for different regions of the spectrum and construct interpolation
functions for them.
We can plot all the spectra included in
sun_elevation.spct for 34 different sun elevation
angles.
We can convert the source_spct object into a
source_mspct to make easier the extraction of individual
spectra.
sun_elevation.mspct <- subset2mspct(sun_elevation.spct, drop.idx = FALSE)Below we plot the spectrum for a sun elevation angle of 45 degrees.
We can also plot them as photon irradiances.
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.
ggplot(sun_elevation.mspct[["tuv-45"]]) +
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()
# each spectrum contains three columns with spectral data that we separate
sun_elevation_tot.mspct <- source_mspct()
sun_elevation_diff.mspct <- source_mspct()
sun_elevation_dir.mspct <- source_mspct()
for (n in names(sun_elevation.mspct)) {
sun_elevation_tot.mspct[[n]] <-
sun_elevation.mspct[[n]][ , c("w.length", "s.e.irrad")]
sun_elevation_diff.mspct[[n]] <-
with(sun_elevation.mspct[[n]],
source_spct(w.length = w.length,
s.e.irrad = s.e.irrad.diff.down))
sun_elevation_dir.mspct[[n]] <-
with(sun_elevation.mspct[[n]],
source_spct(w.length = w.length,
s.e.irrad = s.e.irrad.dir))
}Compute irradiances and the diffuse fraction for different SZAs.
wavebands <- c(list(PAR = PAR()), Plant_bands("sensory"))
wbands_tot.tb <-
q_irrad(sun_elevation_tot.mspct, wavebands, scale.factor = 1e6)
wbands_dir.tb <-
q_irrad(sun_elevation_dir.mspct, wavebands, scale.factor = 1e6)
wbands_diff.tb <-
q_irrad(sun_elevation_diff.mspct, wavebands, scale.factor = 1e6)
wbands.tb <-
data.frame(
SZA = as.numeric(gsub("tuv-", "", wbands_tot.tb$spct.idx)),
SEA = 90 - as.numeric(gsub("tuv-", "", 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"next, we plot of PAR and UV photon irradiances, normalized to sun at the zenith, against solar elevation angle. The proportion of UV-B and to a lesser extent UV-A2 radiation is smaller when the sun is lower above the horizon.
ggplot(wbands.tb, aes(x = SEA)) +
geom_line(aes(y = PAR.tot / max(PAR.tot)), colour = "orange") +
geom_line(aes(y = UVA1.tot / max(UVA1.tot)), colour = "violet") +
geom_line(aes(y = UVA2.tot / max(UVA2.tot)), colour = "purple") +
geom_line(aes(y = UVB.tot / max(UVB.tot)), colour = "black") +
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 we desire to estimate new values in-between observed ones.
# 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]])
}Using PAR as example, we plot the spline as a line and the diffuse fraction values from the TUV model as orange points. (Near an elevation angle of 0 degrees I run TUV simulations at closer values of sun elevation to ensure that the fitted interpolation spline does not overshoot a diffuse fraction of 1 (100%).)
ggplot(wbands.tb, aes(SEA, PAR.diff_fr)) +
geom_hline(yintercept = c(0, 1), linetype = "dashed") +
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)) +
scale_y_continuous(name = expression("Diffuse fraction, "*Q[s] / Q[t]~~(""/1)),
sec.axis = sec_axis(name = expression("Direct fraction, "*Q[d] / Q[t]~~(""/1)),
transform = function(x) {1 - x} )) +
expand_limits(y = 0)
We have created a named list of function definitions, each of them implementing a different interpolation spline.
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.
ggplot(wbands.tb, aes(SEA)) +
geom_hline(yintercept = c(0, 1), linetype = "dashed") +
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)) +
scale_y_continuous(name = expression("Diffuse fraction, "*Q[s] / Q[t]~~(""/1)),
sec.axis = sec_axis(name = expression("Direct fraction, "*Q[d] / Q[t]~~(""/1)),
transform = function(x) {1 - x} )) +
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 computation used here is based on an empirical formula that deviates from simple geometry.
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)
Additional simulated spectral data
Package ‘photobiologyInOut’
supports the import of irradiance and spectral irradiance data from
various file formats as well as direct calls to the Quick TUV
calculator. Data similar to those in sun_elevation.spct
used in the section above, for other conditions or variation in other
variables can be, thus, easily obtained if needed.
References
Aphalo, Pedro J. (2023) High frequency weather data for Viikki, Helsinki, Finland. OSF, .
Madronich, Sasha (n.d.) Tropospheric Ultraviolet and Visible (TUV) Radiation Model. https://www2.acom.ucar.edu/modeling/.
Madronich, Sasha, and Siri Flocke (1999), The role of solar radiation in atmospheric chemistry, in Handbook of Environmental Chemistry, edited by P. Boule, pp. 1-26, Springer-Verlag, Heidelberg.