Skip to contents

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] 525567

The 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 MB

Given 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 TRUE

Function 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 TRUE

Function 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 rows

Function 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.027546

The 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.

autoplot(sun_elevation.spct, 
         annotations = c("-", "peaks"))

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.

autoplot(sun_elevation.mspct[["tuv-45"]], 
         annotations = c("-", "peaks"))

We can also plot them as photon irradiances.

autoplot(sun_elevation.mspct[["tuv-45"]], 
         annotations = c("-", "peaks"), 
         unit.out = "photon")

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.